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Abstract 

Multistationarity in biological systems is a mechanism of cellular decision making. 
In particular, signaling pathways regulated by protein phosphorylation display features 
that facilitate a variety of responses to different biological inputs. The features that 
lead to multistationarity are of particular interest to determine as well as the stability 
properties of the steady states. 

In this paper we determine conditions for the emergence of multistationarity in 
small motifs without feedback that repeatedly occur in signaling pathways. We derive 
an explicit mathematical relationship ip between the concentration of a chemical species 
at steady state and a conserved quantity of the system such as the total amount of 
substrate available. We show that ip determines the number of steady states and 
provides a necessary condition for a steady state to be stable, that is, to be biologically 
attainable. Further, we identify characteristics of the motifs that lead to multistationa- 
rity, and extend the view that multistationarity in signaling pathways arises from 
multisite phosphorylation. 

Our approach relies on mass-action kinetics and the conclusions are drawn in full 
generality without resorting to simulations or random generation of parameters. The 
approach is extensible to other systems. 

Keywords: steady state; kinase; stability; cross-talk; phosphorylation 

1 Introduction 

Multistationarity (the existence of more than one steady state under particular biological 
conditions) in cellular systems can be seen as a mechanism for cellular decision mak- 
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ing. How it arises is therefore fundamental to the understanding of cell signaling, that 
is, the communication of signals to regulate cellular activities and responses. Generally, 
cell signaling involves post-translational modifications of proteins, such as phosphorylation, 
acetylation, or methylation. These modifications change the state of a protein in a discrete 
manner, for example from an active to an inactive state. 

In eukaryotes, reverse phosphorylation is the most frequent form of protein modification 
affecting ~30% of all proteins in humans |T]. Kinases catalyze the transfer of phosphate 
groups to target proteins and phosphatases catalyze the reverse operation. After the com- 
pletion of the human genome project, genome analysis estimated the number of kinases 
to ~ 500 [2j, while the number of phosphatases is smaller by two thirds pQ. Two protein 
phosphatases, PP-1 and PP-2A, account for the vast majority of all phosphatase activity 
[3] with more than 50 PP-1 targets being characterized [3]. 

As a consequence, there is a substantial complexity in the interplay between enzymes 
(kinases and phosphatases) and substrates, exemplified by systems where protein substrates 
use the same catalyzing enzymes (enzyme sharing), and systems where different enzymes 
catalyze the same reaction (enzyme competition). Competition and sharing are general 
examples of cross-talk between motifs. 

The aim of this work is to determine characteristics that lead to multistationarity. 
Following different modeling strategies it has already been shown that feedback in signaling 
networks as well as multisite phosphorylation can both account for multistationarity (5j El 

a. 

We present a mathematical approach for analyzing the steady states of small systems. 
Our method leads to explicit conditions for when multistationarity occurs in terms of rate 
constants and conserved total amounts of substrates and enzymes. Further the approach 
provides means to study the stability of steady states. 

First we present the motifs that we analyze and then we develop the method to deter- 
mine multistationarity and to study stability. The paper concludes with some perspectives 
and discussion. 

2 Motifs 

2.1 Description 

We analyze the motifs shown in Fig. [5j The motifs are referred to as Motif (a)-(l) and 
provide simple abstract representations of known cellular systems. Some examples moti- 
vating our choice of motifs are given in Table [TJ A rich source of examples is found in the 
well-studied MAPK cascades. 

To understand how multistationarity relates to enzyme usage we base our investigation 
on a motif that does not show multistationarity itself. Therefore we build the motifs from 
a one-site phosphorylation cycle which is monostable |16[ [T71 [TBI [T9] and shown in Motif 
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(a) . A specific kinase (phosphatase) catalyzes phosphorylation (dephosphorylation) and 
all modifications can be reversed. In general, protein phosphoforms are denoted by S and 
P, Fig. [5] If one phosphoform is converted into another, an arrow is drawn and the enzyme 
{E or F) catalyzing the reaction is indicated. 

Motifs (a)-(d) cover different possibilities for a one-site modification process. In Motif 

(b) , the same enzyme catalyzes phosphorylation and dephosphorylation. Motifs (c) and 
(d) account for competition between kinases and/or phosphatases to catalyze the same 
modification(s). 

In eukaryotes, phosphorylation of most proteins takes place in more than one site |20| , 
potentially with different biological effects [21] • Combination of two one-site cycles into 
a two-site sequential cycle yields three motifs: (e) all enzymes are different, (f) only one 
kinase but two phosphatases, and (g) one kinase and one phosphatase. By symmetry, 



Motifs Biological phenomena 

(b) A kinase acting also as phosphatase on the same substrate, e.g. HPrK/P 
kinase- phosphatase in Gram positive bacteria [8]. 

(c) ,(d) Several kinases and/or phosphatases acting on the same substrate, e.g. (i) 

Several kinases phosphorylate the alpha subunit eukaryotic initiation factor 
(eIF2a) at Ser51 [9j; (ii) The phosphatases MPK-1 and PTP-SL both modify 
ERK1 

(e) Multi-site phosphorylation by different kinases and phosphatases at each site, 
e.g. (i) Primed kinases, e.g. GSK-3 [TT]; (ii) Aktl is (de)activated through 
three-site sequential (de) phosphorylation by three different kinases (phos- 
phatases) [3]. 

(f) ,(g) Multi-site phosphorylation with the same kinase and/or phosphatase respon- 

sible for all modifications, e.g. (i) Two-site phosphorylation of ERK catalyzed 
by MEK; (ii) Dephosphorylation of ERK2 catalyzed by DUSP6 [12] . 
(h),(i) The same enzyme catalyzing the modification of two different substrates, e.g. 

the kinases ERK1, ERK2 and the kinase products of the p38 pathway catalyze 
phosphorylation of two substrates (the mitogen- and stress-activated protein 
kinase (MSK) 1/2 and the MAP kinase signal-integrating kinase (MNK) 1/2) 

(j),(k),(l) Cascades with several modification steps and substrates, e.g. (i) MAPK 
cascades; (ii) Protein kinase A (PKA) phosphorylates phosphorylase kinase, 
which in turn phosphorylates glycogen phosphorylase (with dephosphorylation 
carried out by the same phosphatase, PP-1, in the two different layers) 
Fig. 7.17], [15]. 



Table 1: Cellular systems represented by Motifs (a)-(l) 
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Motif (f) represents as well a motif with one phosphatase but two kinases. We assume 
for simplicity that both phosphorylation and dephosphorylation proceed in a sequential 
and distributive manner [22j; that is, one site is (de)phosphorylated at a time in a specific 
order. 

Motif (h) represents one-site modification of two substrates that share the same kinase 
but use different phosphatases. This motif represents by symmetry also a system with 
shared phosphatase. If both the kinase and the phosphatase are shared, we obtain Motif 

Finally, two one-site modification cycles can be combined in a cascade motif, where the 
activated substrate of the first cycle acts as the kinase of the next. The interplay between 
enzymes is represented by three cascades: (j) dephosphorylation at each layer uses different 
phosphatases, (k) the phosphatase is not layer specific, and (1) the kinase of the first layer 
catalyzes the modification in the second layer as well. 

2.2 Mathematical modeling 

We assume that any modification S —> S* follows the classical Michaelis-Menten mecha- 
nism in which an intermediate complex Z is formed reversibly but dissociates into product 
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Figure 1: Motifs composed of one or two one-site cycles. Motifs with purple label, and only these, 
admit multiple biologically meaningful steady states. Si and Pj are substrates with i = 0,1,2 
phosphorylated sites. E, Ex, E2 denote kinases, and F, F^Fz phosphatases. In Motif (b) the kinase 
and the phosphatase are the same enzyme. 
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and enzyme G irreversibly: 

S + G ^=± Z — ^ S* + G 

b 

The phosphate donor, generally ATP, is assumed to be in large constant concentration 
and hence embedded into the rate constants. Imposing mass-action kinetics, the species 
concentrations over time can be modeled by a system of polynomial differential equations. 
For example, in Motif (a) the equations are (here E also refers to the concentration of the 
kinase E, and similarly for the other species): 

E=(b E + c E )X - a E ES X = -(b E + c E )X + a E ES 

F = (b F + c F )Y - a F FS 1 Y = -(b F + c F )Y + a F FS l 

S = b E X + c F Y - a E ES Si = c E X + b F Y - a F FS u 

where X (Y) is the intermediate complex formed by the enzyme E (F) and the substrate 
So (Si), and x denotes differentiation of x = x(t) with respect to time. For all motifs, 
there are conservation laws which define time-conserved quantities (total amounts), e.g. 
E + X = and so E = E + X is conserved. The total amounts are fixed by the initial 
concentrations and determine the state space of the dynamical system. Motif (a) has three 
conserved total amounts, namely F = F + Y and S = So + Si + X + Y in addition to that 
of E. 

The steady states of the system are solutions (potentially negative) to the polynomial 
equations obtained by setting all derivatives to zero with the constraints imposed by the 
conservation laws, once total amounts have been fixed. These laws imply that some steady 
state equations are redundant, e.g. either E = or X = can be disregarded. We focus 
on the biologically meaningful steady states (BMSSs), that is, the steady states for which 
all concentrations are non- negative (positive or zero). If at least two BMSSs exist for fixed 
total amounts, then the system is said to be multistationary. 

The specific form of the chemical reactions for Motifs (a)-(l) together with the cor- 
responding systems of differential equations are described in the Supplementary Material 
(SM), attached at the end of the main text. 

3 The Steady State Function cp 

In this section we outline the procedure used to analyze the motifs. Details of the mathe- 
matical analysis are in the SM. 

The system of equations describing the steady states can be reduced substantially by 
elimination of variables [23]. For the motifs considered here, elimination of variables 
implies that the steady states are characterized by a relation S = <p(Y) between the 
concentration of one of the species, typically an intermediate complex Y, and the total 
amount of a substrate S. The concentrations of the other species are given in terms of Y, 
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usually as ratios of polynomials in Y. By imposing all concentrations to be non-negative, 
Y is restricted to a set T of possible values. Further, for any S > 0, there is at least one 
BMSS, that is, S = <p(Y) for some Y in T. The function cp is continuous and differentiable 
in r and depends on the rate constants and the total amounts, except for S. 

The number of BMSSs can be found from the analysis of ip. If (p is strictly increasing 
or decreasing in T, ip is one-to-one and hence, for a given total amount S, there is a 
corresponding unique Y at steady state. Consequently, multistationarity cannot occur 
(Fig. 

Figs. [Bfj and [6J3 show situations where multistationarity occurs. If <p has increasing and 
decreasing parts, or if T is not connected, then Y\ / Y2 with <p(Yi) = (p(Yz) = S might 
exist. Hence, there are at least two BMSSs with the same S. 

These two figures represent substantially different switch responses. In Fig. [6}:, there is 
only one BMSS for low S. An increase of S to S max causes the system to switch to a 'high' 
steady state (high Y) under the assumption that the green steady states in the figure are 
stable. If S is decreased again to S m i n , then the system switches back to a 'low' steady 
state. In Fig. ^p, there is one BMSS for low S. An increase of S keeps the system in the 
first branch of (p and thus it will behave as a monostationary system. 

Interestingly, the derivative (p'(Y) of tp(Y) provides means to determine whether some 
steady states are unstable. Unstable steady states are unattainable under biological con- 
ditions. Specifically, we find that either the regions in which ip is increasing or those in 
which it is decreasing must correspond to unstable steady states, see Section [5] 

In summary, the function <p determines whether multiple BMSSs exist and encodes 
information about the stability of steady states. In the next section we analyze <p for Motifs 
(e) and (f). We show how enzyme sharing in a two-site cycle (f) leads to multistationarity, 



s s s 




(a) (b) (c) 

Figure 2: Possible shapes of <p in T (colored regions: magenta=unstable BMSSs; green= (possible) 
stable BMSSs). (a) ip is increasing and for any s, there is one BMSS (y) such that ip(y) = s. 
(b) T consists of two disconnected regions. For s < S m i n there is one BMSS; for s = S m i n there 
are precisely two; and for s > S m - m there are three; ip is also defined in the white region but some 
concentrations become negative, (c) tp is in part decreasing, in part increasing. For < s < S ma ^ 
there are three BMSSs; for s — S m i n or s = S max , there are two; and for s < S m i n or s > S^nax, 
there is one. 
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as opposed to a two-site cycle with different enzymes (e). A detailed analysis of all motifs 
is given in the SM. 



4 Mono- versus Multistationarity 

4.1 Monostationarity 

Motifs (a)-(e), (h), and (j) have exactly one BMSS for any choice of rate constants and 
total amounts. In all cases, the function ip is increasing in T. The procedure is very similar 
in all cases and is thus only illustrated for Motif (e). We take some effort in explaining the 
details as the procedure might have general applicability. 

Motif (e) consists of three phosphoforms of the substrate, So,Si,S 2 , with subscript 
indicating the number of phosphorylated sites. The chemical reactions of the system are: 

50 + Ei < > Xi — ■>■ Si + Ei Si + E2 < ? X2 — >• S2 + E2 

01, E b2,E 

a l,F Cl f ° 2 > F C 2 F 

51 + Fi Yi — '->■ S + Fi S 2 + F 2 ^Y 2 —^Si + F 2 

We denote the inverse of the Michaelis-Menten constants of Ei by Ki t E = a>i,E/(bi,E + c i,e) 
and of Fi by K^p = a^F/ip^F + (H,f)- The ratio of the catalytic constants of phosphatase 
and kinase is denoted by [i; L = (H,F/ci,E- 

The system has five conserved total amounts, which are assumed to be positive: Four 
for enzymes, Ei = Ei + Xi and Fi = Fi + Yi (i = 1,2), and one for substrate, S = 
So + Si + S2 + Xi + X2 + Y\ + Y%. The steady state equations can be rewritten as 

X 1 = k 1 , e E 1 S Y 1 = K l , F F l S l Xi=A*i*i (1) 

X2 = k 2 ,eE2Si Y 2 = k 2 ,fF2S 2 X 2 = fi 2 Y 2 . 

The last column gives Xi in terms of Yi. The total amounts Ei, Fi give Ei, Fi in terms of 
Yi as well: E t = E t - mYi, Fi = F { - Y { . Further, if E u Fi > then Ei = or F { = 
cannot be solutions to Eq. ( |26[ ). It follows that the concentrations Ei, Fi are positive if and 
only if Yi is in Ti = [0, &) with & = min(Fj, Ei/m). 



We further isolate So, Si from the first row in Eq. (26) and S 2 from the second and 
obtain ^ ^ 

° ki :E (Ei - niYi) ' * K i:F (Fi-Yi) *~ ^ 

for i = 1, 2. Then, S'o, /Si (resp. S2) are non-negative increasing continuous functions of Y\ 
in Ti (resp. I2 in F 2 ). The remaining equation, X2 = k 2: eE 2 Si, gives I2 in terms of Yi: 



(3) 
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The function / is non- negative increasing and continuous in F\. Further, for Y 2 to be in 
T2, it is required that Y\ is in V = [0,^) C Ti with £ = min(£i, f~ l (£, 2 ))- 



Finally, using Eq. (29), we find that X 2 and S 2 are increasing functions of Y\ in T. 
Therefore, using the formulas above, all concentrations at steady state are non-negative if 
and only if Y\ is in T. We conclude that the BMSSs of the system satisfy 

S = So + Si + S 2 + Xi + X 2 + Yi + Y 2 = (p(Yi) 

for Y\ in V. Since 99 is a sum of increasing continuous functions in Y\, then so is <p. 
Additionally, ip(0) = and (f{Yi) tends to infinity as Y\ tends to £. Thus, p has the form 
in Fig. [6^, with a unique Y\ for any given S, that is, there is one BMSS. 

4.2 Multistationarity 

We consider a two-site phosphorylation system with one kinase but different phosphatases 
for each phosphoform, as shown in Motif (f). Multistationarity has been observed numer- 
ically in this system [6]. The system reduces to Motif (e) by setting E\ = E 2 and we use 
the notation introduced previously. The conservation laws are the same with the exception 
that there is only one kinase law, E = E + X\ + X 2 . Define £j = min(i ? j, E/m) and 

r* = [0,6). 



The system of equations to be solved is similar to Eq. (26) with E = E{. Thus, we start 
by writing Xi,E,Fi as functions of Y±,Y 2 . Since E,Fi must be positive at any BMSS, we 
require < Yi < Fi and n{Y\ + [x 2 Y 2 < E. For these values we obtain 

_ ^Y, Y % 

•JO — 7=: — TTT 1 Si 



k>i,e(E - fiiYi - ll 2 Y 2 ) Ki tF (Fi-Yi) 

for i = 1,2, which are non- negative increasing continuous functions of Y^. Using X 2 
k 2 ,eFSi, we obtain Y 2 as a non-negative continuous function of Y\ in T\: 

Y 2 = f(Y^- 



fMi{Kl,F(F 1 -Yi) + K2 t EY 1 ) 



This function resembles that in Eq. (29) except for the quadratic term in the numerator, 
which is a consequence of the conservation law for E involving both Y\ and Y 2 . Further, / 
might not be increasing for all Y\. 

Let T = {Y\ G Ti, such that f(Yi) S T 2 }. Using the formulas derived above, all con- 
centrations at steady state are non-negative if and only if Y\ is in V. Hence, for any 
BMSS, 

S = S + S 1 + S 2 + X 1 +X 2 + Y 1 + Y 2 = y{Y x ) 

with Yi in T. The function tp is continuous with <p(Q) = but T might not be a connected 
interval. 
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Figure 3: The function (p for Motif (f) for different values of F 2 and fixed E > ^\Fi and A > 0. 
Let N — (E — /iiFx)//i2- The values M and a depend on E, Fi and the rate constants (see text). 
The set T is disconnected only when N < F 2 < /(«)■ For large F 2 , <p approaches the black line. 
The vertical bars mark the boundary of F 

Define A = (1 + H2,eI ' k\^f)ijliF\ — E. If A < 0, then / is an increasing function in 
r and we conclude that there is exactly one BMSS. If A > 0, then / has a unique local 
maximum for some a in T\ and all cases in Fig. [6] can occur. By varying the value of F2 
while keeping the other constants fixed, we obtain (Fig. [3]): 

(i) F2 < (E — \x\F\)ln2 (orange): T = [0,ai) with f(a±) = FY The function /, and 
thus cp, are increasing and there is one BMSS (Fig. [6^1). 

(ii) (E — fiiFi)/n2 < F2 < f(a) (green): T = [0, «i) U (0:2,^1) with at < a < 02 and 
f(ai) = f{o.2) = F2. Hence, / is increasing in [0, a{), decreasing in (02,^1) an d 
multistationarity occurs (Fig. [BJd). 

When f(a) < F2 there is an M such that: 

(hi) f(a) < F2 < M (purple): V = [0,^i). The function ip has a decreasing part and 
multistationarity occurs (Fig. [6b) . 

(iv) M < F2 (blue): T = [0, ^1). The function <p is increasing and there is one BMSS 
(Fig. 

4.3 Understanding Multistationarity 

Motifs (f), (g), (i), (k) and (1) exhibit multistationarity for some choices of total amounts 
and rate constants (Fig. [4]). The regions for which multistationarity occurs are detailed in 
the SM. In Motifs (i), (k) and (1) it appears only as in Fig. [6b, while in Motifs (f) and (g) 
both the forms in Fig. [Hja and Fig. [Hj^ occur. 
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It is remarkable that in Motifs (f), (k) and (1) multistationarity occurs for any set of 
rate constants and depends only on the initial conditions (that is, the total amounts). 
Thus, multistationarity can occur in these systems independently of the specific kinetics. 
In contrast, multistationarity in Motif (i) depends on the rate constants and hence not all 
kinetics exhibit multistationarity. The same appears to be the case for Motif (g) |244 [25] . 

The common characteristic of these motifs is that a single enzyme is responsible for 
catalyzing two different substrate modifications, which at the same time are linked (Fig. [4]). 
Indeed, in Motifs (f) and (g) the substrates are linked through Si, which is a modified as 
well as an unmodified substrate for the shared enzyme E. For the Motifs (k) and (1), the 
link is given by Si which is a modified substrate and a kinase, and the common enzymes 
are F and E, respectively. In Motif (i) the kinase E is common and the phosphatase F 
provides the link (or vice versa). In contrast, in Motif (h) an enzyme is responsible for 
two different modifications, but there is no link between the two substrates. Consequently, 
multistationarity cannot be observed. 

Multistationarity can arise from two opposing dynamics acting on the same substrate 



Two-site modification. Shared enzyme: E; Link: Si 




Two substrates. Shared enzyme: E; Link: F 



So 
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HiF > E > fi 2 F 

"Pi l i 2 F>E>n 1 F 



(Hi - A»2)(mi<5i»?2 - £12<52!?l) < 

Two-layer cascade. Shared enzyme: E; Link: Si 




F» Mi.Fi, F 2 large 



Figure 4: Conditions for multistationarity are given. The shared enzyme is marked with a colored 
square; the link is marked with a colored circle; predominant modifications are marked in bold. 
The symbol 3> is short for 'much larger'. 
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(Fig. [4J). For example, in Motif (f), if F\ is much bigger than E and F2 < M, then there 
are multiple BMSSs. Thus, since the amount of phosphatase in the first cycle is much 
larger than the amount of kinase, Si is pushed towards the unmodified form So, while in 
the second cycle S± is driven towards the fully modified form S2 (because F2 < M). 

In Motif (i), provided the conditions on the parameters are fulfilled (Fig. [4|, multista- 
tionarity occurs if either [i\F > E > ^F or H2F < E < [i\F. It implies that in one cycle 
the phosphatase 'wins', while in the other the kinase does. 

5 Stability Analysis 

BMSSs are defined as steady states for which all concentrations are non-negative. However, 
a steady state is biologically attainable only if it is (asymptotically) stable, that is, nearby 
trajectories are attracted to it. We show here for our motifs that if f'(Y) < for some 
steady state (p(Y) = S, then it is unstable. 

5.1 The Jacobian and variable elimination 

For a system of ordinary differential equations in M m , a steady state z is asymptotically sta- 
ble if all eigenvalues of the Jacobian evaluated at z have negative real parts [26} Thm. 1.1.1]. 
Since the Jacobian is a real matrix, the complex eigenvalues come in pairs of conjugates 
and their product is a positive number. If m is odd and all eigenvalues have negative real 
parts, their product, and hence the determinant of the Jacobian, must be negative. If m 
is even and z stable, then the product of the eigenvalues must be positive. Thus, the sign 
of the determinant of the Jacobian provides a necessary condition for a steady state to be 
stable and a sufficient condition for it to be unstable. 

For x = (xi, . . . , x n ) let = (x\, . . . , £), . . . , x n ) (x with Xj removed). We make 
the following observation (see SM for a proof): Let / = (/1, . . . , f n ) : f2 C M" — > R n be a 
differentiable function defined on an open set f2 and z such that f(z) = 0. Assume that Xj 
can be eliminated from the equation /, = in a neighborhood around z; that is, there exists 
a differentiable function tp: ftC?) C IR™" 1 -> M, in VL^ = {x^\x G O}, such that Xj = 
1>(x®) if fi(x) = 0. Define /": -> M n_1 by f k (x) = f k ( Xl , . . . , Xj -i, ip{x), Xj , x n _i) 
for all k ^ i and let J denote the associated Jacobian. Then, the determinant of the 
Jacobian of / at z satisfies 

(-l) l+ ^(z) det(J(^))) = det(J(z)). (4) 

5.2 Unstable steady states 



The relation between the sign of the determinant of the Jacobian and stability, together 
with Eq. (73), lead to a criterion to detect unstable steady states. For each motif, let 
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x = (xi, . . . ,x n ) be the species concentrations, Xi = hi(x) the differential equations and 
Ai = gi(x), . . . , A c = g c (x) the equations for the total amounts. We choose the order of 



the species such that xi 



1, . . . , c, can be isolated from A{ = Qi{x) and the steady state 



equation X{ = becomes redundant. For fixed total amounts, Ax, ... , A c , the steady states 
are the solutions to the system f(x) = of n equations in n variables with fi(x) = gi{x)—Ai 
for i = 1, . . . , c and fi(x) = hi{x) for i = c+ 1, . . . , n. 

Let J(z) denote the Jacobian of / at z. In the SM we prove: If z is a steady state, 
that is, f(z) = 0, and either (i) n — c is even and det(J(z)) < or (ii) n — c is odd and 
det(J(z)) > 0, then z is unstable. The proof relies on the observation made about the 



eigenvalues and Eq. ( 73 ) . 



The function ip of our motifs is derived through successive elimination of variables 
precisely from the system of equations f(x) = 0. Using Eq. (73), the sign of det(J(z)) at 
a steady state z can be traced back from the sign of the derivative of ip (the Jacobian of a 
system with one equation) by considering the equation number (i), the equation variable 
(j), and the sign of dfi/dxj after each elimination. 

To exemplify the procedure, consider Motif (f), where n = 10 and c = 4. The system 
is (see SM for details): 



fl(x 
/ 2 (x 

fs(x 



E + X l +X 2 -E f 6 (x 

Fi + Yl- Fi f 7 (x 

F 2 + Y 2 -F 2 f 8 (x 

S + S 1 + S 2 + X 1 +X 2 + Y 1 + Y 2 -S f 9 (x 

Xi - ki jE ES fio(x 



X 2 — k 2 ^eESi 
X l - nxY x 
X 2 - fi 2 Y 2 
Y 2 - k 2jF F 2 S 2 
Yi - ki pFxSi 



with species x = (E, Fx, F 2 , So, X\, X 2 , Si, S 2 , Y 2 , Yi). The function cp is in terms of 
Y 2 after successive eliminations. Let = ±1 depending on whether the sign of the 
determinant of the Jacobian changes (— ) or not (+) after the fc-th elimination. Then, 
sign((p' (Y 2 )) = (Y\ k efc) sign(det( J(z))), where z is the steady state with z% = Y 2 . 



The order and sign of the eliminations are shown in Table [2j We find that Y\ k et = 1, 
implying that the sign of <p'(Y 2 ) agrees with the sign of the determinant of the Jacobian 
of / evaluated at the corresponding steady state. Since n — c = 6 is even, we conclude 
that the values of Y 2 for which ip is decreasing, that is, p'(Y 2 ) < 0, correspond to unstable 
steady states. Further, it follows that unstable points come between other steady states 
that presumably are stable. 



5.3 Stability in monostationarity motifs 

The Routh-Hurwicz criterion [27] gives sufficient and necessary conditions for the Jacobian 
to have all eigenvalues with negative real parts. Thus, the (asymptotic) stability of a steady 
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(/io,n) 


(2,1,-) 
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5 


(fs,X 2 ) 


(4,2,+) 


+ 











a a) indicates that i,j are the indices of the equation and variable iteratively being 
eliminated and a corresponds to whether fi (fi after substitution of the previous eliminations) is 
increasing (a = +) or decreasing (a = — ) as a function of Xj. 
b obtained as obtained as a(— l) l+J . 

Table 2: Elimination of variables for Motif (f). After each elimination the system / is 
rewritten to correctly determine the sign of dfi/dxj before the next elimination 

state can be determined by this criterion. For the Motifs (a)-(e) and (h) the criterion is 
fulfilled and the unique BMSS is asymptotically stable. We have not been able to determine 
this for Motif (j). 

6 Discussion 

We have investigated small motifs without feedback that account for cross-talk, enzyme 
competition, sharing and specificity in post-translational modification systems and deter- 
mined some features that lead to multistationarity in signaling pathways. 

Bistability, and generally multistability, in biological systems is seen as a mechanism 
of cellular decision making. Compared to systems with a single steady state, the presence 
of multiple stable steady states provides a possible switch between different responses and 
increased robustness with respect to environmental noise. Our study has been driven by 
the observation that biological systems deviate from a one-to-one correspondence between 
enzymes and the modifications they catalyze. This phenomenon, known as cross-talk and 
enzyme sharing, can cause multistationarity and hence be essential for regulating signaling 
systems. 

Our work extends the view of multistationarity as arising from multisite phosphoryla- 
tion [7] to the view that multistationarity is driven by a single enzyme that catalyzes linked 
substrates. Two opposing dynamics acting on the same substrate is a recurrent charac- 
teristic of multistationarity. These observations await a precise mathematical formulation 
and an investigation of its generality. 
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Our approach is conceptually simple and reduces to the study of analytical properties 
of a function ip that relates a conserved total amount and the concentration of a species 
at steady state. The graph (ip(Y),Y) can be seen as a bifurcation diagram with one 
parameter, S. When monostationarity occurs, analysis of tp is quite straightforward, while 
a more in-depth analysis is required when multistationarity occurs. An advantage of this 
approach is that unstable steady states are readily detected from the form of ip. 

The existence of ip is not guaranteed in general. For larger systems, a detailed study 
independently of the rate constants cannot be pursued. However, we have shown that for 
cascades of arbitrary length (extentions of Motif (j)), the function <p exists and properties 
of the full cascade can be derived from properties of the building block, the one-site cycle 
[23]. We are currently working on extending this approach to other systems. 

There are some mathematical characterizations of monostationarity in chemical reac- 
tion networks that relate to our work. The theory of monotone systems [28j characterizes 
systems in which there is only one BMSS and at the same time gives conditions for global 
stability of the BMSS. However, a condition for the theory to be applicable here is that no 
species take part in more than two reactions [29]. This condition is only fulfilled for Motif 
(a) and hence the theory cannot be applied here. 

The only general theory of applicability is that of injective systems |30[ [3T| 132] . The 
motifs that do not allow multiple steady states are in fact injective in the sense of [30| 
when modeled as a continuous flow stirred tank reactor. This fact implies that at most one 
steady state exists [33]. The advantage of this theory is that monostationarity is derived for 
more general kinetics than mass-action [33]. However, when restricted to mass-action, our 
approach is as simple as checking for injectivity and provides in addition simple rational 
functions that enable further comprehensive studies of variation in species concentrations, 
such as stimulus- response curves and signal amplification |23l [35] . 
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Supplementary Material 
A Introduction 

In this Supplementary Material (SM) we provide details of the analysis of multistationarity 
given in the main text, as well as proofs of the results mentioned there. We go through 
all motifs and the reader will easily see that many arguments are repeated as the motifs 
share common structures. We have tried to keep the analysis of each motif independent 
of the analyses of the other motifs. However, the motifs progress from simple one-site 
modification cycles to more complex motifs and the line of thought transpires most easily 
from the simple motifs. It is therefore advisable to study the simple motifs before the more 
complicated motifs. To keep the SM self-contained, Figures 1 and 2 of the main text are 
reproduced here. 

Motif (g) (a futile cycle with two parts) has been studied extensively in the literature. 
Motif (j) (a linear cascade with two layers) was studied for arbitrary length in [23], where 
we showed that the system admits only one steady state. It is briefly covered here for 
completeness. 
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Figure 5: Motifs covered in this work. Figure 1 of the main text. 
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Figure 6: Possible shapes of ip in T (colored regions: magenta=unstable BMSSs; green= (possible) 
stable BMSSs). (a) (p is increasing and for any s, there is one BMSS (y) such that ip(y) = s. 
(b) T consists of two disconnected regions. For s < S* m i n there is one BMSS; for s — S m i n there 
are precisely two; and for s > S m - m there are three; ip is also defined in the white region but some 
concentrations become negative, (c) ip is in part decreasing, in part increasing. For S m i n < s < S max 
there are three BMSSs; for s — S min or s — S max , there are two; and for s < S m - ln or s > 5 max , 
there is one. Figure 2 in the main text. 



A.l Notation and preliminaries 

A continuous and differentiable function with continuous derivative is said to be C . We 
denote by M+ the set of positive real numbers and by M+ the set of non-negative real 
numbers. A rational function is a quotient of two polynomials. An increasing (decreasing) 
function / fulfills f(x) < f(y) (f(x) > f(y)) for x < y, i.e. we take increasing (decreasing) 
to mean strictly increasing (decreasing). The notation x G A means that x belongs to the 
set A. 

We make frequent use of the Implicit Function Theorem in two dimensions to relate 
two variables to each other and to find derivatives of implicit functions. Details about the 
Implicit Function Theorem can be found in text books on functional analysis. 

For x = x(t) a real function of t, we denote by x the derivative of x with respect to t, 
dx/dt. 

A. 2 Constants and variables 

We consider the rate constants to be fixed and positive, i.e. in M + . The constants are 
at,K, c* and those derived from these: At*,*, the inverse of the Michaelis-Menten constants 
and //*, the ratios of the catalytic constants of phosphatase and kinase. Here, to ease the 
presentation, rj* denotes K* t E for a kinase E and (5* denotes K* t p for a phosphatase F. 

The total amounts are likewise considered fixed and positive. Species concentrations 
are considered variables of the system. For example, in X — tjESq, t\ is a fixed unknown 
constant and X, E and So are variables that depend on each other, e.g. X can be considered 
a function of E and So, X = X(E,Sq). That is, the differential equations describing the 
system is a set of polynomials with coefficients in the ring M[a*, fe*, c*, 5*, P*, E*, F*\ or 
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M[r/*, <5*, /x*, S 1 *, P*, E*, i 7 *] (where some constants might not be present in all systems, 
e.g. P is not in Motif (a)) and variables being a finite list of species concentrations, 
E* , F* , 5 1 * , P* , X* , (where some variables might not be present in all systems, e.g. Pq 
is not in Motif (a)). 



A. 3 BMSSs and total amounts 

We only consider the system at steady state, that is all differential equations are put to zero 
and solved for the species concentrations under the constraints imposed by the conservation 
laws. Only solutions to the system of equations for which all species concentrations are 
non-negative, i.e. positive or zero, are of interest. These solutions are called Biologically 
Meaningful Steady States, or BMSSs. Also, we assume that the total amounts are positive 
(non-zero), i.e. enzymes and substrates are always present in the system in some form 
(e.g. in some phosphoform or intermediate complex). 

One can prove for each specific motif that if all total amounts are positive, then all 
species concentrations at a BMSS are positive as well (i.e. if a BMSS exists then all species 
concentrations must be positive). It follows that if 1) the species concentrations are non- 
negative (i.e. positive or zero), 2) the total amounts are positive, and 3) the conservation 
laws and the steady state equations are fulfilled, then the species concentrations constitute 
in fact a BMSS and hence are positive. This observation is frequently used in the following. 



A.4 Method 

For all motifs we follow the same procedure. We take the set of differential equations 
describing the systems together with the equations for the total amounts (the conservation 
laws) and solve for one variable, say an intermediate complex Y. Specifically, we choose 
one equation for a total amount, e.g. S = So + S\ + X + Y, and use the other equations 
(differential equations and equations for total amounts) to provide expressions for the 
variables at steady state in terms of Y, i.e. we find expressions such that Sq = Sq(Y), 
Si = Si(Y) and X = X(Y) are functions of Y. These functions only depend on Y, the 
rate constants and the total amounts, excluding S. The functions are substituted into the 
expression for S to obtain a function ip(Y) that relates Y to S, i.e. S = <p(Y). The analytic 
form of <p determines how many BMSSs the system has for a given set of rate constants 
and total amounts. For example, if ip is increasing there can only be one Y corresponding 
to a given S. See Figure [6] for illustration. 

We allow Y to be zero in which case S = f(Y) also is zero (it follows from the con- 
struction of (p for each motif). This is done for simplicity as we otherwise would have to 



consider the limit of <p(Y) for Y — > in each case (see Section A. 3 about positive and 
non- negative concentrations). Some variables (depending on the motif) are not allowed to 
be zero as this could lead to division by zero. 
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A. 5 Technicalities 

Some manipulations of the equations are left to the reader if they only involve standard 
elimination, rewriting and differentiation techniques. In general, it is a good idea to do 
the calculations yourself (either by hand or in Mathematica , Maple or similar) as many 
details are left out due to space and readability constraints. We have done all derivations 
in Mathematica™ and checked them manually. Proofs of propositions are collected in an 
appendix to keep the presentation simpler. Equation numbers are to the immediate right 
of the equations and only equations that are used later are numbered. 



B Steady states of the motifs 
B.l One-site phosphorylation cycles 

Motif (a). We recently used this motif as the building block of linear cascades and 
showed that it admits only one BMSS [23J. We follow the approach taken in |23j and 
summarize it below. Essentially this is the approach taken for all motifs. 
The chemical reactions of the motif are 

S + E ^=± X S 1 +E S x +F zr^~ Y S + F 
b E b F 

The corresponding system of differential equations is: 

S = b E X + c F Y - a E ESo S t = c E X + b F Y - a F FS l (5) 

E = (b E + c E )X — a E ESo X = -{b E + c E )X + a E ES (6) 

F = (b F + c F )Y-a F FS 1 Y = -{b F + c F )Y + a F FS l . (7) 

It follows that E + X = 0, F + Y = and S + Si +X + Y = 0. Hence, the conservation 
laws are given by 

E = E + X, F = F + Y, S = So + Si+X + Y. 

Due to the conservation laws some equations are redundant, for example F = and Y = 
are the same equation. If fixed total amounts are given, we have to solve a system of 6 
equations in 6 variables consisting of the equations for the total amounts and, for instance, 
equations ^5J)-([7p. From the latter equations, we obtain the equivalent system 

X - t]ESo = 0, Y — bFS x = 0, X — /jY = (8) 

with constants rj = a E / (b E + c E ), 5 = a F / (b F + c F ) and \i = c F /c E . 

It follows that at steady state E = E — [iY and F = F — Y . Note that if both E, F are 
positive (i.e. non-zero), then we cannot have E = or F = at steady state (for example 
if E = 0, then X = from @, hence E = 0). Let f = rsAn{F,E/n) and V = [0,0- The 
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variables E, F are both positive and Y is non-negative only if Y is in T. Further it follows 
from Q that 

So = — ==- > and Si = — = . 

n(E-fiY) 8{F-Y) 

(Note that division by zero does not occur as E = E — /uY and F = F — Y are both greater 
than zero.) These two functions are non-negative, increasing and C 1 for Y € T. When Y 
tends to £, one of them tends to +00, while the other remains bounded. 

Let A = fiF - E, so that T = [0,F) if A < 0, and V = [0,E/(j) if A > 0. Using the 
conservation law for S we obtain: 

Result 1 (Motif (a)). Let a one-site modification cycle be given with positive total amounts 
S, E, F. Then the system has a unique BMSS. Specifically, the BMSS satisfies S = f{Y) 
for Y in T = [0, £), £ = mm(F,E/fi), where 

S = tp(Y) = - Y + -J- + (1 + u)Y 

^ V ; n(E - (iY) 5(F - Y) K W 

is an increasing rational C 1 function which tends to infinity as Y tends to £ and fulfills 
<p(0) = 0. 

Remark 1. Since ip is a rational function, the equation S = (f(Y) can be written in 
polynomial form by elimination of denominators. In the present case 

p(Y) = fi5Y(F - Y) + r]Y(E - fiY) + 7]6{(1 + fi)Y -S)(F - Y)(E - nY), 

which is a third degree polynomial in Y. Note that p(0) < 0, > 0, p(() < with 
< £ < C = max(i ? , E/fj,), and p(Y) tends to +00 as Y tends to +00; hence p(Y) has 
three positive roots. However, only the first root is in T, and it corresponds to the only 
BMSS of the system. In some systems, several positive roots are BMSSs. This remark is 
applicable in all cases below whenever ip is a rational function. 

Remark 2. We argued above that E and F are non-zero if E and F are positive. This 
ensures that we do not divide by zero when constructing ip. For the remaining variables 



we only need to ensure non- negativity, cf. Section A. 3 For the other motifs, we make use 
of similar reasoning. 

Motif (b). Consider the system in Motif (b) where the two catalyzing enzymes are the 
same. The chemical reactions of the system are given by 

S + E X Si + E S! + E Y S + E 
The corresponding system of differential equations is: 
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S = biX + c 2 Y - ai ES Si = dX + b 2 Y - a 2 ES! (9) 

E = (b 1 + c 1 )X + (b 2 + c 2 )Y-a 1 ES -a 2 ES 1 X = -(h + c x )X + ai ES (10) 

Y = -(b 2 + c 2 )Y + a 2 ES v (11) 

We find that E + X + Y = and So + Si + X + Y = 0. Hence, the conservation laws 
are given by 

E = E + X + Y, S = S + Si + X + Y. 

If the total amounts are given, we need to solve a system of 5 equations in 5 variables 
consisting of those for the total amounts and, for instance, equations ([£])-( 11). From the 
latter equations we obtain an equivalent system given by 

X - rjESo = 0, Y — 5ESi = 0, X — fiY = 

with constants rj = a%/(bi + ci), 5 = a 2 /(b 2 + c 2 ) and fj, = c 2 /c\. Note that if E > 0, then 
E = cannot be a solution of the steady state equations and thus we require E > at any 
BMSS. 

It follows that E = E — (fi + 1)Y. For E > and Y > we have that Y must be in 
r = [0,15/(1 + /i)). Further, we find 

<->! 



^-( M + 1)Y)' ^-(/x + l)Y) 

These functions are continuous and increasing for Y in V. In addition, all concentrations 
So, S\,X, E, F are non- negative as functions of Y if and only if Y G V. Using the conser- 
vation law for S we obtain: 

Result 2 (Motif (b)). Let a one-site modification cycle be given with one enzyme acting 
as kinase as well as phosphatase. Further, assume that the total amounts S, E are positive. 
Then, the system has a unique BMSS. 

Specifically, the BMSS satisfies S = <p(Y) for Y inT = [0, £) ; £ = E/(l + fj,), where 

S = ip(Y) = _ ^ Y + -= + (1 + u)Y 

^ ' n(E-(p + l)Y) 5{E — (fi + 1) Y) 1 N 

is an increasing, rational C 1 function which tends to infinity as Y tends to £ and fulfills 
^(0) = 0. 



Motif (c). Consider now a one-site modification cycle with two competing kinases. The 
chemical reactions of the system are given by (k = 1, 2): 
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S + E k X k — Si + E k S l + F Y S + F 
The corresponding system of differential equations is: 

E k = (bf + 4 )X k ~ a%E k S X k = -{4 + c$)X k + a* E k S (12) 

F = (b F + c F )Y - a F FS 1 Y = -(b F +c F )Y + a F FS 1 (13) 

S = bfX 1 +b F X 2 + c F Y -a F E 1 S Q -a F E 2 S Q Si = cf X 1 + cf X 2 + b F Y - a F FS! (14) 

with k = 1,2. It follows that X k + E k = 0, Y + F = and S + Si + X\ + X 2 + Y = 0, 
leading to the conserved total amounts (k = 1,2): 

E fc = £ fc + X fc , F = F + Y, S = S + Si + X x + X 2 + Y. 

Therefore, if total amounts are given, the steady states of the system are solutions to a 
system of 8 equations in 8 variables consisting of the equations for the total amounts (4 
equations) and, for instance, equations (12) for k = 1,2, (13) and (14). From the latter 
equations we obtain an equivalent system given by 

X k -rj k E k S = 0, Y-5FSi = 0, cfX 1 +c%X 2 -c F Y = (15) 

for k = 1,2, with constants r} k = a k /(b^ + c^) and 5 = a F /(b F + c F ). Note that if 
E k ,F > 0, then neither E k = nor F = are solutions of the steady state equations. 
Thus, E k ,F > at any BMSS. 

From the total amounts we have that E k = E k — X k and F = F — Y . Hence, for 
E, F > and Y > we must have that < X k < E k , < Y < F at BMSS. We find 

Si = — r, (16) 
S(F-Y) V 7 

which is increasing in Y. Likewise, we obtain the following relation (with non-zero deno- 
minators) 

X l _ a _ X 2 

TUfa-Xi) ~ °~ V2(E~2 -X 2 )' 

and it follows that 

Xl = x (x 2 ) = _ VlElX2 . 

rj 2 (E 2 - X 2 ) + r\\X 2 

This function is increasing and C 1 for X 2 in T 2 = [0, E 2 ). Additionally, < X\ = 4>x(X 2 ) < 
E± for X 2 in T 2 . Note that So is also expressed as a non- negative, increasing C 1 function 

of x 2 e r 2 . 
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From the last equation in (15) we obtain 



Y = MX2) = (cf/c F W x (X 2 ) + (cf /c F )X 2 

which is a non-negative, increasing, rational C 1 function of X 2 for X 2 G IV 

If F is in (/>y(T2), then the condition Y = 4>y{X 2 ) < F sets a more restrictive upper 
bound to X 2 than E 2 does. The functions </>v,0y ar e defined for X 2 = E 2 , and we have 
y (£ 2 ) = (cf/c^)^! + (c F /c F )E 2 . Let a = (cf/c^l + (c F /c F )E 2 , f = min(F,a), 
f = ^(f) (which is well-defined) and T = [0,0 C T 2 . Then for X 2 G I\ we have Y < F. 



Since Si is a non-negative, increasing C function of Y (for Y < F, (16)), it is also a 
non-negative, increasing function of X 2 G T. In conclusion, all steady state concentrations 
are non- negative if and only if X 2 G T. Additionally, either So or Si tend to infinity as X 2 
approaches £. Using the conservation law for S we obtain: 

Result 3 (Motif (c)). Let a one-site modification cycle be given with two different competing 
kinases and one phosphatase. Further, assume that positive total amounts S , E±, E 2 , F are 
given. Then the system has a unique BMSS. 

Specifically, the BMSS satisfies S = (f{X 2 ) for X 2 inT = [0, £), where 

S = <p(X 2 ) = X2 + ,J y{ ^ y ^ + 4x(X 2 ) + X 2 + MX2) 
n 2 (E 2 -X 2 ) 5{F - 4>y{X 2 )) 

is an increasing, rational C 1 function which tends to infinity as X 2 tends to £ and fulfills 

■m - 0. 

Motif (d). Consider a one-site modification cycle with two competing kinases and two 
competing phosphatases. The chemical reactions of the system are given by 



c 



E af r F 



S + E 1 ^^X 1 ^+S 1 +E 1 S 1 +F 1 ^^Y 1 ^+S + F 1 

So + E 2 < t X 2 > S\ + E 2 S\ + F 2 < t Y 2 >■ So + F 2 

The corresponding system of differential equations is: 

X k = -(b% + c^)X k + a^E k S (17) E k = (6jf + cf )X fe - af £ fc 5 

n = — C&fc + cf )Y fc + 4 F^i (18) F fc - (&£ + cf )Y fe - af F fc 5x 

S = 6f Xi + cf Y - afEiSo + b%X 2 + c%Y 2 - a%E 2 S 
^ = cf X x + 6f Yi - of F1S1 + cf X 2 + b%Y 2 - af F 2 5i (19) 
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with k = 1,2. It follows that X k + E k = 0, Y k + F k = and S + S x +X x +X 2 +Y x +Y 2 = 0, 
leading to the following conserved total amounts (k = 1, 2): 

Therefore, if total amounts are given, we have to solve a system of 10 equations in 10 
variables consisting of the equations for the total amounts (5 equations) and, for instance, 



equations (17), (18) for k = 1, 2 and ( 19 ). From the latter equations we obtain an equivalent 
system given by 

X k - Vk E k S o = 0, Y k -6 k F k S x = 0, cf X x + cf X 2 - cf Y 1 - cf Y F = (20) 

for k = 1, 2 with constants r/ k = af / (&f + cf ) and 5 k = af /(&f + cf). 

If E k ,F k > 0, then E k ,F k ^ at steady state. As for the previous motif, we have 
E k = E k — X k and F k = F k — Y k ; hence < X k < E k and < Y k < F k is required for 
any BMSS. The concentration So can be expressed in two different ways as increasing C 1 
functions: As a function of X x and as a function of X 2 . When these two expressions are 
equated we obtain the relation (similar to the relation obtained for the previous motif) 

mE x x 2 



X 1 = <t >x (X 2 



rj 2 (E 2 - X 2 ) + 771X2 



It is a non- negative, increasing C 1 function for X 2 G [0, E 2 ) such that E x = E\ — 4>x{X 2 ) 
also is a positive function of X 2 G [O,^)- Note that (ftx(E 2 ) = E\ is well-defined. 

Similarly, S\ can be expressed as increasing C 1 functions of Y\ and of Y 2 , respectively, 
which provide the relation 

Yi = 4>y{Y 2 ) 



8 2 {F 2 -Y 2 ) + 5 1 Y 2 

It is a non-negative, increasing C 1 function for Y 2 G [0,^2) and F\ = F\ — 4>y(Y 2 ) is a 
positive function of Y 2 G [0,^2). 

Finally, the last relation in ( |20[ ) gives 

cf 0x(a 2 ) + cf x 2 = cf^ y (y 2 ) + cfy 2 . 

The left side is an increasing C 1 function in X 2 , the right side an increasing C 1 function 
in Y 2 which tends to infinity as Y 2 tends to infinity. Hence, there exists an increasing C 1 
function g(X 2 ) = Y 2 defined on X 2 G [0,^2) relating X 2 to Y 2 . 

In summary, the concentrations E X ,E 2 , X\ , Sq are non- negative functions of X 2 if and 
only if X 2 G [O,-^)- The concentrations F x , F 2 , Y x , S\ are non-negative if and only if 
Y 2 G [0,-F2)- Hence, to ensure that all concentrations are non-negative, we require y2 = 
g(X 2 ) < F 2 . Since g is increasing, it is bounded above by g{E 2 ) (well-defined). Therefore, 
let r = [0,E 2 ) if g(E 2 ) < F 2 and T = [0, g^ 1 (F 2 )) otherwise. Using the conservation law 
for S we obtain: 
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Result 4 (Motif (d)). Let a one-site modification cycle with two competing kinases and two 
competing phosphatases be given. Further, assume that the total amounts S,Ei,E2,F\,F2 
are positive. Then, the system has a unique BMSS. 

Specifically, the BMSS satisfies S = f(X 2 ) for X2 inT = [0, £), where 

S = <p(X 2 ) = X2 + - ^ X f + 4>x(X 2 ) + X 2 + 4^{jg(X 2 )) + g(X 2 ) 
ri2{E 2 - X 2 ) d2{F2g{X 2 )) 

is an increasing C 1 function which tends to infinity as X2 tends to £ and fulfills ty?(0) = 0. 

The function g is not rational, hence neither ip is rational. 

B.2 Two-site phosphorylation cycles 

Motif (e). First we consider a two-site phosphorylation system in which modifications are 
carried out by different kinases and phosphatases for each phosphoform. For simplicity, we 
assume that both phosphorylation and dephosphorylation occur sequentially. The chemical 
reactions of the system are: 

ci^ qE 

So + Ei < > X\ >• S% + Ei Si + E2 < ? X2 >■ S 2 + E2 

bf b$ 



j 2 

2 f cf a 2 



Si + Fi ~ — > Yi — U- S + Fi S 2 + F 2 ~ — " Y 2 — S 1 + F 2 
6f bf 

The differential equations describing the system are: 

X k = a%E k S k -i - (fef + cf )X k E k = (&f + cf )X k - af S fc 5 fc _i (21) 

Y k = af F k S k - (&f + cf )F fe , F fc = (&f + cf )Y fc - af F k S k (22) 

S = &f*i + cf Fi - of £i5 S 2 - cf X 2 + &f Y 2 - af F 2 S 2 (23) 

Si = cfX 1 + b%X 2 + b(Y t + cf Y 2 - (af Fi + a%E 2 )S 1 (24) 

with k = 1,2. We have X fe + £ fc = 0, Y fc + F fc = and 5 + 5i + 5 2 + Xi+X 2 + Yi +^2 = 0, 
which lead to the following conserved total amounts (k = 1,2): 

E k = E k + X k , F k = F k + Y k , S = So + Si + S 2 + Xi + X 2 + Y 1 + Y 2 . (25) 

Therefore, if total amounts are given, we have to solve a system of 11 equations in 11 



variables consisting of those in (25) and, for instance, equations (21)-(24). From the latter 
equations, we obtain 

X k = r) k E k S k -i, Y k = 5 k F k S k , k = 1,2 (26) 
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with constants r\k = a^/(b^ + cj?) and S k = a k /(b^ + c k ). Equation (22) and (23) for 
k = 2 give C2X2 — C2Y2 = 0. This relation together with equation (24), (21 ) for k = 2, and 
(22) for k = 1 give cfX\ — c±Y\ = 0. Therefore, we have that 

X k = /j, k Y k , fj, k = c k /c k . 

Let A k = fj, k F k - E k , 6c = Tam.(F k ,E k /iik) and T k = [0,£ fe )- Note that if E k ,F k > 0, 
then at steady state E k , F k 7^ 0. Since E k = E k — fJ, k Y k and F k = F k — Y k at steady state, 
any BMSS must satisfy Y k E T k to ensure non- negativity of X k , Y k and positivity of E k , Fk 
as functions of Y k for each k. 

For fixed values of S2, X2, Y2, the steady state values of the remaining variables satisfy 
the steady state equations of Motif (a) (a one-site phosphorylation cycle) with species 
5o, Si, Xi, Yi, Ex, Fi and total amounts Ei, F\ and S — S2 — X2 — Y2. Therefore, using 
Result [TJ the BMSSs of the system satisfy 

S = ipi(Y 1 ) + (S 2 + X2 + Y 2 ) 

with Yi G Ti. Here <pi denotes the function <p in Result [Tj 

Using the second equality in ( 26 ) for k = 2 together with the conservation law for F2 , 
we obtain 

MY2) = S 2 + X 2 + Y 2 = - - Y2 - + 112Y2 + Y 2 , (27) 

d 2 {F 2 - Y 2 ) 

which is a non-negative, increasing C 1 function of Y 2 G T2. Consequently, the BMSSs of 
the system satisfy the relation 



S = ( P1 {Y 1 ) + ip 2 (Y 2 ) 



(28) 



for Yi G Ti and Y 2 ET 2 . The right hand side is an increasing function in both variables. 

From the equation X2 = rj 2 E 2 Si together with Si expressed as a function of Yi (similar 
to that of S2 in equation (27)) we obtain 



X2 = m(E2-x 2 )Si 



m(E 2 -x 2 )Yi 

6i(F 1 -Y 1 ) 



Rewriting this equation yields 
m E 2Yi 



X 2 



5i(Fi-Yi)+ V2 Yi 



and hence Yz = f(Yi) 



m E 2Y 1 



fj, 2 Si(Fi — Yi) + ii2m Y \ 



(29) 



The function / is an increasing function, which is non-negative and C 1 for Y\ G Yi. Addi- 
tionally, / is defined for Y\ = F\ with f(F{) = E 2 j ' ji 2 . 
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By substitution of f(Yi) into ( [28| ), the BMSSs of the system satisfy 

S = <p{Y x ) = <p x (Y 1 ) + <p 2 tf(y 1 )) 

for Y\ E T\ and f(Yi) £ T 2 . Since / is increasing, continuous and /(0) = 0, this condition is 
equivalent to Y\ € Y = [0, £) with £ = min(£i, Z" 1 ^)) (note that / _1 (£2) is well-defined). 
The function </? is a rational function, which is increasing and C 1 in T and either (f\ or 
(f2 ° f tends to infinity as Y\ tends to £. Additionally, the BMSS concentrations of all 
other species derived from the formulas above are non-negative functions of Y\ if and only 
if Y\ G T. 

Using the conservation law for S we obtain: 

Result 5 (Motif (e)). Let a two-site phosphorylation cycle be given with different kinases 
and phosphatases. Further assume that the total amounts E k , F k , S, k = 1, 2 are positive. 
Then, the system has a unique BMSS. 

Specifically, the BMSS satisfies S = (p(Y\) for Y\ inT = [0,£) ; where 

S = <p(Y 1 ) = <p 1 (Y 1 ) + ip 2 (f(Y 1 )) 

is an increasing rational C 1 function, which tends to infinity as Y\ tends to £ and fulfills 
99(0) = 0. 



Motif (f). Next, we consider a two-site phosphorylation system where phosphorylation is 
catalyzed by the same kinase at both sites but dephosphorylation is catalyzed by different 
phosphatases. Again, we assume sequential (de)phosphorylation. The chemical reactions 
of the system are: 

5 + E X 1 -^-^ Si+E Si+E X 2 -^-^ S 2 + E 

51 + Fi t~~^" Y x 5 +F X S 2 + F 2 Y 2 S x + F 2 
The differential equations describing the system are the following: 

5 = bfX x + c F Y x - afESo S 2 = cf X 2 + b F Y 2 - a F F 2 S 2 (31) 
Y k = -(&£ + c k) Y k + 4 F k S k X k = -(6f + cf )X k + af £S fe _i (32) 
E = (bf + cf)X 1 + (b F + c F )X 2 -(a F So+a F S 1 )E F k = (b F + c F )Y k - a F k ' F k S k (33) 

51 = cfX 1 + b F X 2 + b F Y 1 + c F 2 Y 2 - (of F 1 + a F E)S 1 (30) 

with k = 1,2. The conservation laws are given by 

E = E + X x + X 2 , F k = F k + Y k , S = S + S 1 + S 2 + X 1 + X 2 + Y 1 +Y 2 , 
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k = 1,2. Therefore, if total amounts are given, the system to be solved consists of 10 
equations in 10 variables which are the equations for the total amounts and, for instance, 



equations (30)- (33) for k = 1, 2. From (32) and (33) we obtain 



X k = rjkESk-x, Y k = 5 k F k S k , k = l,2 

with constants rj k = af /(&f + cf) and 5 k = a k /(b k + c^). Proceeding as in the previous 
system, we find that 

X k = /ZfeYfc, fi k = c k /c k . 

Let A fc = fi k F k -E, £ k = mm(F k ,E/fj, k ) and Tfc = [0,^). If E,F k > 0, then at steady 
state E, F k ^ 0. From E = E — X\ — X2 and F k = F k — Y k we see that positivity of E, F k 
and non-negativity of Y k requires that (at least) Y k £ T k . If Y k £ T k , then X k is also a 
non- negative, increasing function of Y k . 

The situation resembles the situation of the previous system where catalysis is mediated 
by two different kinases. In both systems we have X2 = 11-2X2 and S 2 = j=r an d so 

02(-T2— 12) 

the concentration of 5*2 is a function of Y2. It follows that 

<P2(Y 2 ) = S 2 + X 2 + Y 2 = - - - + ii 2 Y 2 + Y 2 

02(^2 - Y2) 

is a non- negative, increasing continuous function of Y 2 G T2- 

For a fixed value of Y 2 , the steady state values of the remaining variables (except X2 
and S2) satisfy the steady state equations of a one-site phosphorylation cycle with species 
So, Si, Xi, Yi, E, F\ and total amounts E — X2 = E — fJL 2 Y 2 , F\ and S — (p 2 (Y 2 ). Using 
Result [T] with (pi(-,Y 2 ) denoting ip for a fixed Y 2 , the BMSSs satisfy the relation 

S = Vl{Y x ,Y 2 ) + ? 2 {Y 2 ) 

for any Y 2 £ T2 and < Y\ < min(Fi,(E — fJ, 2 Y 2 )/[ii). Under these conditions, all 
concentrations of the other chemical species are non- negative functions of Y±,Y 2 . Note 
that the total amount of enzyme E — H2Y2 is part of the function ip\ and hence (pi depends 
on Y 2 . Indeed, we have 

DO = = • 

Tji(E - \i\Y\ - \x 2 Y-i) 

The function (pi is increasing in Y\ and in Y 2 . 

The equation X2 = rj 2 ESi combined with S± expressed as a function of Y\ provide 
(after isolation of X 2 ) the following relation at steady state 

Y 2 = m) =x 2 /, 2 - ^-^ yi 



H 2 Si(Ft -Y{) + ^mYi 
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(A) (B) (C) 



Figure 7: Different behaviors of the function f(Y x ) = Yi- The region T x is marked in 
green. (A) corresponds to Proposition [l] (i) with Y x G T x = [0,F X ). (B) and (C) differ in 
whether (B) Y x e Y x = [0,E/m) or (C) Y x € Ti = %F X )\ in both cases /(B/mi) = 0. If 
^2 > /(«) (the top point of the curve) then Proposition [I] (h)(c) applies. If F2 < f(a) 
then Proposition [l] (h)(b) applies for (B). In (C), if F2 < f(F x ) (the dashed line) then 
Proposition 1 (ii) (a) applies, whereas Proposition 1 (ii)(b) applies if f(F x ) < F2 < f(a). 



This function resembles that in (29) except from the quadratic term in the numerator 



which gives / a very different analytic form from that in ( 29 ) . 

The function / is C 1 and takes non-negative values for Y x < £1. Therefore, a BMSS 
satisfies 

S = ip{Y x ) = MYi, f(Yt)) + ¥> 2 (/(n)), (34) 

for Y\ S Ti, such that /(ii) G T2- Since both ipi,tp2 are increasing functions, the behavior 
of </? needs to be understood from the behavior of /. 
If we let 

t = {y 1 € ri|/(ii) g r 2 }, 

then the concentrations of all the chemical species are non-negative when expressed in 
term of Y x , if and only if Y x G T. If Y x = 0, then /(0) = and it follows that 93(0) = 0; 
in particular, it follows that G T. Note that ip is C 1 in T; however T might not be a 
connected interval as we will see below. 

The function / is C 1 and non- negative for Y x EF X . Depending on whether it is monotone 
or not, different forms of ip are expected. These behaviors can be found by computing the 
derivative of / (see Figure [7]) . Let 

A=(l + m /6 x )pL X F x -E. 

Observe that if A < then A x < 0. 

Proposition 1. The following statements hold: 

(i) If A < 0, then f is an increasing C 1 function ofY x 6 T x = [0, F x ). Further, we have 
r = [0, min(i ? i, /~ 1 (£ 2 )) (with / _1 (£ 2 ) set to +00 if not defined). 
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(ii) If A > 0, then there is a £ T\ = [0, £1) such that f'(a) = 0, / is increasing in [0, a) 
and decreasing in (a,£i). /n i/tis case, we obtain: 

(a) If E — jUi-Fi — ^2-^2 > 0, i/ien i/iere is ai < a sttca £aai / is increasing and C 1 
inT = [0, Qi) and /(ai) = F2. 

(ft) If E — [iiFi — H2F2 < and F2 < f(a), then F = [0, a±) U («2, £1) ai < a < 
«2 and /(«i) = f(ct 2 ) = F2. Hence, f is increasing in [0,ai) and decreasing in 

("2,£l)- 

fcj IfF 2 >f(a), then V =[0,^). 

All possibilities are covered in the proposition. Indeed, the condition in (ii) (c) implies 
E — [j,iFi — H2F2 < 0: If not, according to (ii)(a) there would be ot\ < a such that 
F2 = /(cki) < /( a ) < F2, which is a contradiction. 

In the cases (i) and (ii)(a), / is an increasing C 1 function in a connected interval 



r = [0,0- Hence, by composition of functions, ip in (34) is also an increasing C 1 function 
of Y± G r with ip(0) = 0. Additionally, in both cases (p tends to infinity as Y\ tends to ^ 
since either I2 = tends to F2 or Yi to F\. We conclude that there is exactly one 

BMSS in each of these cases. 

The cases (ii) (b) and (ii) (c) require further analysis. In case (ii)(b), let 

T' = [0, ai ), r" = (a 2 ,&). 

The values a%, 02 correspond to the values of Y\ for which f{Y\) = F2 and hence they are 
the zeros of the denominator of 52: 

q Y 2 Y 2 



S 2 (F 2 -Y 2 ) 6 2 (F 2 -f(Yx)) 

In r', / is increasing and therefore (p is also increasing. It tends to infinity as Y\ tends to 
a%. Hence, for any value of S, there is a BMSS corresponding to a Y\ € V . 

In r", / is a decreasing function and hence it is uncertain when/whether cp is increasing. 
However, we find that (A) when Y\ tends to 0.2 from the right, the function (P2 f tends to 
+00, while Tp x = (p\(-, /(•)) is bounded; (B) when Y\ tends to £1 from the left, the function 
ip2 o / is bounded, while Tp-^ tends to +00 (in fact, either S\ or Sq expressed in terms of Y\ 
does). It follows that in the interval T" = (02, £1), t ne function tp starts decreasing from 
infinity and ends increasing towards infinity. By continuity, there exists a minimum .Smm 
of <p in this interval. We see that for 5 > S m i n , at least two values of Y\ satisfy tp(Y{) = S; 
hence at least two BMSSs exist in this interval. All together, we conclude that at least 
three BMSSs occur in this case, one in T' and two in T" . Note that when S = S m i n then 
there are at least two (not at least three) as the two in T" coincide. 

Finally, let us consider case (ii) (c), that is the case A > 0, F2 > f(ct) and consequently 
r = [0,£i). The function <p is increasing at Y\ = and tends to infinity as Y\ approaches 
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For fixed values of F\ and E with A > 0, multistationarity occurs for certain values of F2 
and S. The figure shows the function / for Ai < colored according to the the shape of ip as 
function of F2 ■ 

_ Y 2 
If F2 > M (blue region) there is only one BMSS. 
If M > F 2 > f(a) (purple region) there is mul- 
tistationarity, (ii)(c) in Proposition 1. If f(a) > 
F2 > — A1//X2 (green region), then multistation- 
arity also occurs, (ii)b. Finally for — A1//12 > F2 
(orange region) there is only one BMSS, (ii)(a). If 
Ai > 0, then the orange region cannot occur. 




Figure 8: Different shapes of tp depending on E,Fi,F%. 



£1. Multistationarity can only occur in the form of Figure |6]^c) and there will be different 
values of Y\ corresponding to the same value of S, only if the function p is not always 
increasing; that is if the derivative of p has more than two zeros. Equivalently, if there 
exists Y\ for which p'{Yi) < 0, then we are guaranteed multistationarity. This result is 
stated in the following proposition: 

Proposition 2. Assume that A > 0. Then there exists a value M > /(a) such that, for 
F2 > M, the derivative of p is positive ( except potentially for a finite number of points 
where it is zero). Hence, ip is increasing. Further for M > F2 > f(ct), there exist values 
ofYi for which <p'(Yi) < 0. 

We conclude that for all F2 > M, there is exactly one BMSS for any value of S, while 
for all M > F2 > f{oe) there is multistationarity for certain values of S. In particular, 
multistationarity occurs for S satisfying S m j n < S < S max , where S m in 1S the smallest local 
minimum (excluding 0) and S max the largest local maximum of p (excluding infinity), cf. 
Figure [6^c). 

Based on numerical examples we have made the following observation: When multi- 
stationarity occurs, the function p has one local minima and one local maxima, f3\ < 02, 
resulting in at most three BMSSs. When F2 increases, /3% increases, while fa decreases. 
For some F2, fa = fa and the decreasing part of p is lost and p becomes increasing. In 
fact, the value of F2 for which fa = fa is M. 

All together, this implies that for fixed values of E,F±, the value of F2 determines 
whether values of S for which multistationarity occurs exist. Figure [8] shows how the 
number of steady states changes with i^- 

We conclude that: 
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Result 6 (Motif (f)). Let a two-site phosphorylation cycle be given with one kinase cata- 
lyzing phosphorylation at both sites and two different phosphatases catalyzing dephospho- 
rylation. Further, assume that positive total amounts E, Fi, F2, S are given. 

Then, the BMSSs satisfy S = <p(Yi) for Y\ 6 T = [0,£), where ip is a C 1 function with 
<p(0) = 0. 

Let A = (1 + r]2/5i)niFi — E. Then we have: 

(i) If A < or E — /-i±Fi — /U2-F2 > then there is a unique BMSS. 

(ii) If A > then there exists values of F2 and S such that the system has more than one 
BMSS. Further, there is an upper bound to F2 for which multistationarity can occur. 
This bound is independent of S. 



Motif (g). Next, we consider a two-site phosphorylation cycle, where phosphorylation is 
catalyzed by the same kinase at both sites and dephosphorylation is catalyzed by the same 
phosphatase at both sites. Again, we assume sequential phosphorylation. The chemical 
reactions of the system are: 



S + E: 



bf 



S! + E 



Si + E: 



b§ 



X 2 



S 2 + E 



Sx + F: 



S + F 



S 2 + F. 



bf 



St + F 



The differential equations describing the system are the following: 



S 2 = cfx 2 + b^Y 2 - a£ FS 2 (36) 



72 = C 2 j\ 2 f U 2 1 2 



S = bf Xi + cf Yi - afESo 

E=(bf + cf )Xi + (bf + cf)X 2 - (afS Q + afS^E X k = -{bf + cf)X k + afES k ^ (37) 

F=(bf + cf)Y 1 + (bf + cf)Y 2 -(afS 1 +afS 2 )F Y k = -(bf + cf)Y k + af FS k (38) 

^ = cfX x + bfX 2 + bfYx + cfY 2 - (afF + afE)S 1 (35) 

with k = 1, 2. As in the previous system, the conservation laws are given by 

E = E + X 1 + X 2 , F = F + Y 1 +Y 2 , S = S + Si + S 2 + X x + X 2 + Y X + Y 2 . 

Therefore, if total amounts are given, the system to be solved consists of 9 equations in 
9 variables which are the equations for the total amounts and, for instance, equations 
(35)|- ([38]). From (37) and (g we obtain 

X k = Vk ES k - 1} Y k = 5 k FS k , k = l,2 (39) 

with constants r] k = af / (bf + cf) and S k = o,f/(bf + cf). From (36) and (38) with k = 2, 



and from (35), (36), (38) with k = 1,2 and (|37j) with k = 2,we obtain 

X k = n k Y k , n k = cf/cf. 
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If E, F > then E, F ^ at steady state. We isolate E, F from the corresponding 
conservation laws and write them as functions of Yi,Y 2 . Since the denominators are non- 
zero, we find 

a = Mjjl a = Y 2 

° Vl (E - ^ - fi 2 Y 2 y 2 S 2 (F-Y 1 -Y 2 y 



From Yi = 5 1 FS 1 = 5 1 (F-Y 1 - Y 2 )S l we obtain Y x 
X 2 = r] 2 ESi = rj 2 (E - niYi - n 2 Y 2 )Si, we find that 



<SiSi(F-y 2 ) 

1+<5iSi • 



Now from the equality 



Y 2 



We can therefore write all concentrations as functions of Si . Let 

Pi(Si) = E + 6 1 (E-fi 1 F)S 1 

p 2 {Si) = n 2 F + m {n 2 F-E)S 1 

Pz{Si) = \i 2 + fj, 2 (Si + r] 2 )Si + 5ir/ 2 (^2 - Hi)Sf. 



Then we have 

5 1 Sip 2 (Si) 



Y 



P3{Sl) 



Yo 



T] 2 SlPl(Sl 



So 



fJ>i5iSip 2 (Si 



S-2 



inSipijSi 

5 2 p 2 {Si) 



This leads to the following relation for the BMSSs 

S = <p(Si) = S + S 1 + S 2 + X 1 + X 2 + Y 1 +Y 2 

with 93(0) = and where each term is considered a function of Si, as defined above. For 

this motif, a 93-function in terms of an intermediate complex cannot be obtained: The 

relationship between Si and Y x or Y 2 is in general not invertible. 

We next seek to determine the values of Si that lead to non-negative concentrations of 

the species; that is to determine the region T, where Si is defined. From (39) we see that 

if Yjfc, So, Si, S 2 are non-negative, then so are E, F, at steady state. The roots of pi,p 2 

are 

E fi 2 F 



6 



Si(niF-EY r\ 2 (E - fi 2 F) 

respectively. Since pi,p 2 are in the denominators of Sq,S 2 , these polynomials are not 
allowed to vanish. Let Ai = /iiF — E and A2 = [i 2 F — E. Then 

(i) If A 2 < 0, then p 2 (5i) > if and only if Si < £ 2 . If A 2 > 0, then p 2 (^i) > for all 
Si > 0. 
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(ii) If Ai < 0, then pi(Si) > for all Si > 0. Otherwise, if Ai > 0, then pi{Si) > if 
and only if Si < £i. 

The polynomial P3 (Si) has degree 2 and is positive whenever Si is positive and H2 > Hi- 
If hi < Mi) then there is exactly one positive root £. Therefore, we have the following 
situations: 

A. If fi2 > Hi (A2 > Ai), we require pi(Si),p2(Si) > 0. Three different scenarios occur: 

1) Ai > and A 2 > 0, i.e. E < m~F < H2F, F = [0,fi). 

2) Ai < and A 2 > 0, i.e. mF <E< /jqF, F = [0, +00). 

3) Ai < and A 2 < 0, i.e. n{F < fi 2 F <E,F = [(),&)■ 

B. If fxi > Hi (^1 > A 2 ), we have 

1) Ai > and A 2 > 0, i.e. E < /x 2 F < HlF, F = [0,fi) for £1 = min(£,£i). 

2) Ai ^0 and A 2 < 0, i.e. H2F < E < HlF, F = [0,0u(f,+oo) for f = min(fi,f 2 ,0 
and £ = max(£i, 

3) Ai < and A 2 < 0, i.e. /i 2 F < HlF <E,F= [0,6) for & = min(f,£ a ). 

Observe that in all cases, the function tp is C 1 and non-negative in F. In the cases A.l, 
A. 3, B.l and B.3, when Si approaches the upper limit of F, then tp tends to infinity, since 
at least one of Y\, Y2, So, <S 2 does. In case A. 2 the function tends to infinity as Si tends to 
infinity. In all these cases, multistationarity would appear in the form of Figure (6^c), that 
is, if should decrease for some values of Si. 

In case B.2, F is not connected so we let F = F' U F" . When Si tends to £, then tp 
tends to infinity. It follows that for every S, there is one BMSS located in F' . Additionally, 
when Si approaches £ from the right, tp also tends to infinity, implying that (p comes down 
from infinity in F" . When Si tends to infinity, tp tends to infinity. This implies that in 
case B.2 the function resembles that in Figure |6^b) , potentially with more increasing and 
decreasing parts. In any case, when S is large, multistationarity occurs and there are at 
least three steady states. 

Let us consider the derivatives with respect to Si of the following summands in tp: 

d(Yi + Y 2 ) = H2{m~E + 81H2F) + 25ir]2Hi{H2 - fn)FSi + 772^1(^2 - ^i)A 2 S 2 
dSi P3 (Si)i 
8(hiYi + H2Y2) _ nlimE + S1H1F) + 25ir]2H2(H2 - Hi)ESi - r] 2 5 2 H2(H2 - A*i)Ai5? 



dSi p 3 (S 



1 2 



dS = 81H1H2E F + 2rj 2 EA 2 Si - / mr ?2 ^A 1 A 2 5 1 2 
dSi (r]iH2 Pi{Si)) 2 

8S 2 = mH2~E F - 2 m 5iH2F/±iSi - rjSiAiA^Sl 
dSi~ (52P2(Si)) 2 
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In case A. 2 these derivatives are all positive since Ai < 0, A2 > and ^2 — A*i > 0. 
Hence, there is exactly one steady state. For the remaining cases, one can always find 
combinations of parameters for which the function (p is decreasing for some value of S\. 
Thus, further analysis of the derivatives or the function ip itself is required. For example, in 
case A. 3, the following choices of numerical values provide multiple steady states: F = 3, 
E = 10, Si = 10, S 2 = 100, m = 1,H2 = 3, m = 0.002, f] 2 = 100. 

We conclude that: 

Result 7 (Motif (g)). Let a two-site phosphorylation cycle be given with phosphorylation 
at both sites being catalyzed by the same kinase and dephosphorylation at both sites being 
catalyzed by the same phosphatase. Further, assume that the total amounts E, F, S are 
positive. Then, any BMSSs satisfy S = <p(Si) for Si € T (with T defined above) where if 
is a C 1 function with y(0) = 0. 
Additionally, 

(i) If niF < E < H2F, then for any total amount S there is a unique BMSS. 

(ii) If H2F < E < niF , then there exists a value S m i n such that for any S > S m i n the 
system has at least three BMSSs. 

B.3 Modification of two different substrates 

Motif (h). In this system, two cycles are connected through a joint catalyzing kinase. 
The chemical reactions of the system are: 

Sq + E ^=± Xi -5-*- Si+E P + E ^=± X 2 Pi + E 

bf bf 



Si + F x T~~^" Yi S + Fi Pi + F 2 Y 2 -i* P + F- 



6f bF 



2 



2 



The differential equations describing the system are the following: 



So 


= bfXi + cf Yi - of ES 




Po 


= b 2 E X 2 + cl-Y 2 -a 2 E EP 




Si 


= b(Yi+ cfXi - af F1S1 


(40) 


Pi 


= biY 2 +ciX 2 -a(F 2 Pi 


(43) 


Xi 


= -(6f +c?)X 1 +afESo 


(41) 


x 2 


= -(6f + cf)X 2 +afEP 


(44) 


Yi 


= -(b( + c()Y 1+ a(FiSi 


(42) 


Y 2 


= ~(b 2 +c^)Y 2 + a^F 2 Pi 


(45) 


Fi 


= (b( + cf )Yi - af fiSx 




F 2 


= (6f + cf)Y 2 -afF 2 P 1 




E 


= {bf + c f)Xi + (bi + ci)X 2 


- (af So + a 2 E P )E. 









The conservation laws are given by 
E = E + Xi+X 2 , F k = F k + Y k , S = S + S 1 + X 1 + Y 1 , P = P + Pi + X 2 + Y 2 , 
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with k = 1,2. Therefore, if total amounts are given, the system to be solved consists of 11 
equations in 11 variables which are the equations for the total amounts (5 equations) and, 



for instance, equations (40)-(45). From (41), (42), (44) and (45) we obtain 



X 1 = rnES Q , X 2 = i l2 EP , Y 1 = SiF 1 Si, Y 2 = S 2 F 2 P h 



with constants = af /(b^ + c%) and 5k = + c%). From these equations, (40 ) and 



(43) we obtain 

X k = HkYk, (J>k = c k/ c k- 

If E, Fk > 0, then at steady state E, F^ ^ 0. Since Fk = F — Y}~ we require < Yfc < F^ 
for any BMSS. Using the conservation laws for P and F 2 we have 

r] 2 E 6 2 (F 2 - Y 2 ) 
and it follows that E is an increasing C 1 function of Y 2 

fi 2 S 2 Y 2 (F 2 -Y 2 ) 



E = g{Y 2 ) 



V2(5 2 P(F 2 - Y 2 ) -Y 2 - 5 2 (fi 2 + 1)Y 2 {F 2 - Y 2 )) 



The numerator is non- negative for < Y 2 < F 2 . The denominator is a degree two polyno- 
mial in Y 2 with positive independent and leading coefficients. For Y 2 = F 2 , the denominator 
is negative and thus has two positive real roots. It follows that for E to be non-negative 
we require Y 2 G [0, £2) where £2 < F 2 is the first positive root of the denominator. We have 
that g(0) = and that g(Y 2 ) goes to infinity as Y 2 tends to £2- To sum up, for < Y 2 < £ 2 , 
the steady state values of F 2 , Po, Pi, X 2 , E are non- negative as well. 
Using the conservation law for E, it follows that 

E = E + X 1 +X 2 = g(Y 2 ) + f , 1 Y 1 + ^Y 2 , (46) 

and Y 2 is a decreasing C 1 function of Y\, h(Y\) = Y 2 , defined on [0,E/fii] such that 
h(E / m) = and h(0) < £ 2 . Consequently, g{h(Y\)) = E is a decreasing C 1 function 
defined on [0, E/^xi] and since E > 0, we require Y\ G [0, E/fii). 

Let T = [0,£i) with £1 = min(Fi, E/fii) and consider the conservation law for S. We 
obtain 

V19KKY)) Ji(Fi-yi) 

Then, ip is an increasing C 1 function defined on T with image M+ and all steady state 
concentrations are non-negative if and only if Y\ G V. Therefore, we have shown: 
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Result 8 (Motif (h)). Let two one-site modification cycles with joint kinase and distinct 
phosphatases be given. Further, assume that the total amounts S,P,E,F k , k = 1,2, are 
positive. Then, the system has a unique BMSS. 

Specifically, the BMSS satisfies S = (p(Yx) for Y\ E T = [0,£i), where 

5 = y(yi)= t^W + tt^ v . +viY l + Y 1 

is an increasing C 1 function which tends to infinity as Y\ tends to £i and fulfills <p{0) = 0. 



Remark 3. In equation (46), if we isolate Y\ instead of Y2, we would obtain a relationship 
between Y± and Y2, Y\ = h(Y2) and the steady state relation in terms of Y2, S = ^(i^)- 
In this case, Y2 6 where £ = /i(£i) and £' is the pre- image of E of the increasing 

function H2Y2 + 9(^2)- If £1 = E/ fjL\, then £ = 0. 

Motif (i). In this system, two cycles are connected through a joint catalyzing kinase and 
a joint catalyzing phosphatase. The chemical reactions of the system are: 

E CL^ E 

S + E ^=± X x Si + E P + E ^=± X 2 Pi + E 

of C F af C F 

Sx+F ~ — *• Y x — ^ S + F Px + F ~ — " Y 2 — ^ P + F 
The differential equations describing the system are the following: 



So 


= bfXi + cf Yx - 


- afES 






Pa 


= bfX 2 - 


f cfY 2 - 


afEPo 




Si 


= b(Y 1 + c?X 1 - 


- afFSi 




(47) 


Pi 


-fefr 2 4 


- cf * 2 - 


afFPi 


(50) 


Xx 


= -(6f +cf)X l 


+ afES 




(48) 


X 2 


= -(bf- 




f af £P 


(51) 


Yx 


= -(b( + c() Yl 


+ af FSi 




(49) 


Y 2 


= -(6f- 


H <£)Y a 4 


- af FPi 


(52) 


F 


= (bf + cf)Y 1 + 


(bf+cf)Y 2 - 


(af Si 4 


a£Px)F 












E 


= (bf + cf)X x 4 


(bf + cf)X 2 


- (af So 


+ afPo)E. 













The conservation laws are given by 

E = E + Xx + X 2 , F = F + Yx + Y 2 , S = S + Sx + Xx + Yx, P = P + Px + X 2 + Y 2 . 

If total amounts are given, the system to be solved consists of 10 equations in 10 variables 
which are chosen to be the equations for the total amounts (4 equations) and equations 



(47)- (52). As usual we derive 

Xx = m ES , X 2 =V2EP , Yx = 6xFSx, Y 2 = 5 2 FPx, X k = fi k Y k , 
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k = 1,2, with constants rjk = af /{bf + cf), 5k = af / (&f + cf) and = cf /cf. 

Let = min(F, £ , //u/ c ), = 1,2. For fixed Y% E [0, £2)) the equations above provide 
the steady state equations corresponding to a one-site cycle with species E, F, So, Si, X%, Y\ 
and total amounts S,E — ^2^2-, F — Y2. Analogously, for fixed Y\ E [0,£i), we obtain the 
equations for a one-site cycle with species E, F, Pq, Pi, X2, Y 2 and total amounts P,E — 
HiYi,F — Y\. Therefore, using Result [I] we obtain that the steady states are solutions to 
the system 

S = MYi,y 2 ) = ^ — + Yx + Ml y 1 + y 1 

j]i[E - H1Y1 - H2Y2) d 1 {F-Y 1 -Y 2 ) 

P = MYi,Y 2 ) = ^ — + Y ' 2 +v*Y 3 + Y 2 

r) 2 {E - mYi - H2Y2) 6 2 {F -Y1-Y2) 

for any Y\,Y2 and for the solution to be a BMSS we require that Y\ + Y2 < F and 
li{Yi + H2Y2 < E. All species concentrations derived with Yi,Yz fulfilling these conditions 
are non-negative. 

Note that (pi, ip2 are increasing functions of both Y\ and Y2. Therefore, using (pi for a 
fixed S, there is a decreasing C 1 function /(I2) = Yi, defined for Y2 E [0, Further, since 
/(I2) is the steady state value of Y\ in the first cycle with total amounts S,E— H2Y2, F—Y2, 
the derived concentrations E, F, So, Si, X\ are also non- negative. Thus, so are Pq,Pi,X2 
fory 2 E [0,6)- 

Let r = [0, £2) so that all concentrations are non- negative if and only if Y2 E T. We 
conclude that the steady states of the system are described by a C 1 function 

P = V(Y 2 ) = (p2(f(Y 2 ),Y 2 ) 

defined for Y2 E T. 

The behavior of the function <p determines the presence/absence of multiple steady 
states. Note that when Y2 tends to £2; then Yi = /(I2) tends to zero and <p tends to 
infinity, implying that (p increases as we approach the upper limit of T. Additionally, at 
Y2 = 0, /(Y2) is finite, <p(0) = and hence the function is increasing at as well (it 
is positive for Y2 > 0). We conclude that the existence of at least one steady state is 
guaranteed and multistationarity occurs if ^'(^2) < f° r some Y2 E T. 

To proceed, let ip = (tpi, ^2) : M 2 — )■ M 2 and let J^u be the Jacobian matrix of ifi, that 

is, 

J1I1 - 

Then, a simple observation (proved in Proposition [7]) shows that ^'(Y^) < if and only if 
det(J^(/(y 2 ), Y 2 )) < 0. Therefore, if det(J^(Yi, Y 2 )) < for some values E,F,Y X ,Y 2 , then 
for total amounts S = (pi(Yi, Y2) and P = 99 2 (Yi, Y 2 ), the system exhibits multistationarity. 
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Proposition 3. Let a = (fix — /i 2 )(^io"in 2 - and T = {(Yx,Y 2 ) € M 2 | Y\ + Y 2 < 

F, Ml y x + /x 2 Y 2 < £}. 

ft) lf°>0, then det(J i ,{Y 1 ,Y 2 )) > /or a// (Yi,r 2 ) € F. 

(zzj Assume that a < 0. // (^aj fii — fJ- 2 > 0, and either fi 2 F > E or E > [i\F ', or (b) 
H2 — /ii > 0, and either fiiF > E or E > fi 2 F, ^en det( J^(YI, Y 2 )) > /or a// 

(Yi,Y 2 )er. 



(raj If a < and either [i\F > E > \i 2 F or [i 2 F > E > fi\F, then there exist values 
(Yx,Y 2 ) G T, suc/i i/iat det(J^,(Yi,F 2 )) < 0. 

The proof of this proposition is found in Appendix |D) Using it, we derive the following 
result: 

Result 9 (Motif (i)). Let two one-site modification cycles with joint kinase and joint 
phosphatase be given. Further, assume that the total amounts S, P, E, F are positive. Then, 
the BMSSs satisfy P = <p(Y 2 ) f or Y 2 *E T, where <~p is a C 1 function which tends to infinity 
as Y 2 tends to £ 2 and fulfills </?(0) = 0. 

Let a = — /i 2 )(/ii<5in 2 - /U 2 <5 2 ni). Then: 

• The function ip is always increasing if either (i) a > or (ii) a < and either (a) 
fJ-i — fJ-2 > 0, together with fi 2 F > E or E > fi±F, or (b) [i 2 — \i\ > 0, together with 
HiF >EorE> ^ 2 F. 



• If a < and either \i\F > E > fi 2 F or fi 2 F > E > n\F , then, there exist values 
Y 2 E r for which <p'(Y2) < 0. Hence multistationarity occurs. In this case the total 
amounts S, P are required to be large. 

B.4 Cascade motifs 

Motif (j). We consider here the combination of two one-site modification cycles in a 
cascade motif with a specific phosphatase acting in each layer. The chemical reactions of 
the system are: 

CL^ qE ^E 

Sq + E < - X\ > S\ + E Pq + Si < > X 2 *- P\ + S\ 

bf bf 

Si + Fx Yx S + Fx Pi + F 2 Y 2 P + F 2 

bf bf 

The differential equations describing the system are the following: 
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S 


= bfx 1 + cf Y 1 


- afESo 


Xi 


= of £S - (bf 4 






(53) 


A 


= bfX 2 + cfY 2 


— af S\Po 


X 2 


= of SiPo - (bf - 


f Cf )*2 




(54) 


E 


= (b? + cf)X 1 - 


- afESo 


Y 


= of - (bf - 


fcf)Y! 




(55) 


A 


= (bf + cf)Y 1 - 


-ofFiSi 


Y 2 


= 4F 2 P 1 -(bi- 


fcf)Y 2 




(56) 


F 2 


= (6f + cf)Y 2 - 


- of F 2 P 1 


A 


= bfY 2 + cfX 2 - 


- of P 2 Pi 




(57) 








Si 


= cfX l + bf Y x 4 


-(6f + cf)X 2 - 


- (af F 1 - 


hafP )Si. (58) 



The conservation laws are given by 
E = E + X 1 , F k = F k + Y k , S = S + S 1 +X 1 + X 2 + Y 1 , 



P = P + Pi + X 2 + Y 2 



for k = 1,2, If total amounts are given, the system to be solved consists of 11 equations 
in 11 variables which are the equations for the total amounts and, for instance, equations 
(53)-(58). As usual we obtain 

X 1 = rjiESo, X 2 = mSiPo, 



Y 1 = S1F1S1, Y 2 = 5 2 F 2 Pi 



with constants r\ k 



+ 4) and S k 



a k/( b k +4)- Equations ph,d57b (k = 2) 



and (54), (55), (58), (k = 1) provide the relations 

X k = fi k Y k , fi k 



c k l c k 



for k = 1,2. 

The concentrations E,F\,F 2 , X\ , X 2 are solved in terms of Y\ , Y 2 from the total amounts 
E, F\,F 2 and the two relations above. If E, F k > 0, then E, F k ^ at steady state. Hence, 
any BMSS satisfies < Y k < F k , < Y x < E/fn. 

Let £i = min(i ? i, E/fii). Then we find 



Ml*l 



and S 



Yi 



miE-vxYt) 6 1 (F 1 -Y 1 y 

which are non- negative, increasing, continuous functions of Y\ E [0,^i). It follows that 

i Pl (Y 1 ) = S + S 1 + X 1 + Y 1 

is a non-negative, increasing, continuous function of Y\ 6 [0,£i). Also, it tends to infinity 
as Y\ tends to £i and thus the image of ipi over [0,£i) is From the conservation law 
for S we find 



Y 2 = /(Yi 



1 

M2 



(S-MYi)), 
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which is a decreasing function in Y\ G [0, £1). The concentration Y2 is non- negative provided 
¥>i(li) < S. Let £2 = l Pi 1 (S) (fi is increasing, hence invertible). Then Y2, Sq, S±, X\ are 
all functions of Y\ and non-negative provided Y\ G [0,^]- We have that f{^i) = and 
/(0) = S / ji2- It follows from the inverse function theorem that / can be inverted so that 
there exists a continuous decreasing function Y\ = 5(^2) for Y2 £ [0, S / ^2] with g(0) = £2 
and (7(577/2) = 0. Then So, Si, Xi, 11 are all functions of I2 and non-negative provided 
Y2 £ [0, S/fi2\- Note that if S > then Si = is not a solution of the steady state 
equations and hence any BMSS satisfies Si > 0. Since Si = g(Y2) / {5\{F\ — (/(I2))) we 
require Y2 < S/^2- 

Let £ = min(F2, S///2) and T = [0,£). Next we find the relations 

P = — ~pr 1 and, Pi - 



mSi 5 2 (F2-Y2) 

which are non- negative continuous functions of Y2 £ T. Since Si is increasing in Y\ and 
thus decreasing in Y2, we have that both Pq, P\ are increasing in Y2. 

We have seen that all concentrations at steady state are non-negative if and only if 
Y2 G r. Further, when I2 tends to £, either Po or -Pi tend to infinity. Using the conservation 
law for P we obtain: 

Result 10 (Motif (j)). Let a cascade of one- site modification cycles be given with positive 
total amounts S, P, E, F\,F2- Then the system has a unique BMSS. Specifically, the BMSS 
satisfies P = ^(12) for Y2 inT = [0, £), where 

is an increasing C 1 function, which tends to infinity as I2 tends to £ and fulfills ip(0) = 0. 

Motif (k). We consider here the combination of two one-site modification cycles in a 
cascade motif where the phosphatase is not layer specific; that is the same phosphatase 
acts in both layers. The chemical reactions of the system are: 



,E 



S + E ^± X l S!+E Po + Si X 2 -i- Pi + S*i 



6f 6f 



Si + F t~~^ li S + F P\+ F ^=± Y 2 P + F 

The differential equations describing the system are the following: 
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Sq — 


bfx 1 


+ cfYi 


- afESo 






= afESo 


~(bf 4 




(60) 


Po = 


bfX 2 


+ cfr 2 


— afS\Po 




x 2 


= af S\Po 


"(6f- 




(61) 


E = 


(bf+ 


cf )Xt - 


- afESo 




Y 


= afFS x 


"(tf + 


cf)Y 


(62) 


F = 


(bf+ 


<$)Yx 4 


-(bi + 4)Y 2 - 


(afS 1 +aiP 1 )F 


Y 2 


= af FPi 




c F 2 )Y 2 


(63) 


Si = 


cfX, 


+ b(Y 1 


+ (bf + d)X 2 


- (of F + afP )S 1 


(59) A 


= bfY 2 + 


cfX 2 - 




(64) 



The conservation laws are given by 
E = E + X 1 , F = F + Y 1 + Y 2 , S = S + S 1 + X 1 + X 2 + Y 1 , P = P + Pi + X 2 + Y 2 . 
If total amounts are given, the system to be solved consists of 10 equations in 10 variables 



which are the equations for the total amounts and, for instance, equations (59)-(64). From 
the latter equations, we obtain 

X 1 = m ES , X 2 = m S 1 P , Y 1 = S 1 FS 1 , Y 2 = 5 2 FP 1 



with constants rjk = af /(frf + cf) and 5k 



and (|59j) , (JSTJ) and @ (jfe = 1) provide 



= cf /cf 



Equations (63), (64) (k = 2) 



for k = 1,2. 

Note that if E, F > 0, then E, F ^ at steady state. Therefore, any BMSS satisfies 
< Yk < F and Y\ < Ej[i\. If S\ = 0, then it follows from the steady state equations 
that So = X\ = X 2 = Y\ = and thus S = 0. Hence, if S > we have that Si > and so 
Y 2 = X 2 /fi 2 < S/n^fov any BMSS. _ _ 

Let £2 = min(5//X2, F), £1 = min(P//ii, F) and = [0, = 1,2. For fixed 

Y 2 £ T 2 , let £1(12) = min(P/^i,P - Y 2 ) and note that £i(Y 2 ) < fr. For F 2 fixed, the 
steady states satisfy the steady state equations of a one-site phosphorylation cycle with 
species So, Si,X\,Y\, E, F and positive total amounts E, F — Y 2 and 5 — ii 2 Y 2 (these are 
independent of P). Using Result [l] with ipi(Yi,Y 2 ) denoting <p(Y\) for the fixed Y 2 , the 
BMSSs satisfy the relation 

S = $(Y 1 ,Y 2 ) = <pi(Yi,Y 2 ) + n 2 Y 2 

for < Y% < £,i(Y 2 ) < £1. Under these conditions, the concentrations So, Si,E, X\,Y\, 
X 2 ,Y 2 , F are non- negative. Note that 



So 



m {E-ix x Y x 



Si 



Yx 

5 1 (F-Y 1 -Y 2 y 
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(with non-zero denominators) and (pi(Yi,Y 2 ) = So + Si + X\ + Y\, such that only Si 
depends on Y2. Since the BMSS is unique in the one-site cycle, it follows that for every 
Y 2 G T2 there exists Y\ G Ti (Y± is non-zero provided S > 0) satisfying S = ^(Yi,!^)- 
Thus, there is a function / so that for any BMSS, 

Yi = f(Y 2 ). (65) 

The function / is C 1 and decreasing for Y2 G T2- Indeed, this follows from the implicit 
function theorem and the fact that <3> is C 1 in T2, the derivative with respect to Y± is 
positive and the derivative with respect to Y 2 is positive too. This can be checked by direct 
computation. 

By construction, we have < /(I2) < £1(^2) for every Y2. Using Remark[TJ the value 
Y\ = f(Y 2 ) is the first positive root of the polynomial in Y\ (for fixed Y 2 ) obtained from 
<J>(Yi, Y2) by elimination of denominators. Also, if Y2 = then Y\ = /(0) is the BMSS value 
of a one-site phosphorylation system with total amounts S, E, F and hence it is positive. 

Consider now the conservation law for P: 

D H2Y2 li2h(F-Y 1 -Y 2 )Y 2 Y 2 Y 2 

rn = — = — , 



mSi V2Y1 5 2 F s 2 (F-Y 1 -Y 2 ) 

and hence P = ^2(^15^2) is a function of Y\ and Y 2 , which is C 1 with respect to each 
variable. Further, each summand Po,Pi,X2 and Y2 of P is non- negative if Y 2 £ T2 and so 
all concentrations at steady state are non- negative if and only if Y 2 £ T2 ■ Accordingly, we 

let r = r 2 . 

The BMSSs are characterized by the relation 

P = v{Y 2 ) = V2{f{Y 2 ),Y 2 ). 

The function <p is C 1 for Y 2 G T. Further, ip(0) = and as Y 2 tends ^2, f(Y 2 ) tends to zero 
and hence <^(Y 2 ) tends to +00. Consequently, there is at least one BMSS. 

Determination of the number of BMSSs follows from the behavior of the function ip in 
r. Since T is a connected interval, multistationarity can only occur if ip'(Jf 2 ) < for some 
value Y 2 G V (Figure [6^c)). In this case multiple BMSSs occur for P m i n < P < P max for 
some -P m in < -Pmax- Analysis of the behavior of this function is not straightforward but 
still some general conclusions can be derived. 

Clearly Y 2 + 112Y2 is an increasing function of Y 2 . The derivatives of Pq, Pi with respect 
to Y 2 are 

dPo li>2&\ , t! v , ov x fivfl? v\\ a 5P i F-f + f'Y 2 

W2=^p U{F ~ f ~ 2Y2) - fY2{F ~ Y2)) and m = ^-/-y 2 f 

< 0. 



Either or |^ must be negative for (p'(Y 2 
Let Ai = mF - E. 
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Proposition 4. Let a = min(5', P) and = min(l + fix, /U2/2). If Ai > and 

a„(F -E/nx)>a 

then §^(^2); §^-(^2) > for all Y2 G T. Consequently, multistationarity cannot occur. 

Finally, let us analyze the situation where S is large. 
Proposition 5. The following statements hold: 

(i) Fix Y2 > F — E I Then, as S tends to +00, ip'Qfo) becomes positive, (p'iY^) > 0. 

(ii) Assume that Ai > 0. Let A 2 = 27 m^l^E^iF - E) -fo^i/^Mi-F- (^281 +772(1 + 
fj,2))E) 3 . If A2 < 0, i/ien c/? decreases for some Y2 < F — E j ' [i\ as S becomes large. 

Implicitly Proposition [5] assumes that P is large too: In (i), Y2 is fixed and S becomes 
large restricting the possible values of Y\ and P. In (ii), the same is in play. Thus, 
contradictory conclusions cannot be reached from the two propositions. For example, if 
F 3> E, then A2 < and Proposition (5^ii) guarantees that for S (and thus P) large, 
multistationarity exists. If F S> E and P, S fixed, Proposition [4] ensures monostationarity. 
Also, it follows from Proposition [5] (i) that if F — E/{i\ < 0, multistationarity for S large 
cannot occur. 

Multistationarity can also occur if Ai < 0. For instance, the total amounts F = 4, E = 
W,S = 50, P = 195 and rate constants /i 2 = 15, /ii = 2,772 = 771 = #i = 1,82 = 0.1 
produce three steady states with Yz = 2.25,2.78,3.26, respectively. However, contrary to 
the situation with the reversed inequality, Ai > 0, multistationarity cannot occur for S 
large (Proposition [5](i)). 

We have the following result: 

Result 11 (Motif (k)). Let a cascade be given with one phosphatase acting in both layers. 
Further, assume that the total amounts E, F, S, P are positive. 

(i) If F ^> E and S large, then there exist values P m \ a < P max such that for all P m m < 
P < -P max the system admits more than one BMSS. 

(ii) IfF — E/fii >~aja il with a = mm(S , P) and = min(l + m, ^2/^), then multista- 
tionarity cannot occur. 

(Hi) If F — E/fj,\ < 0, then multistationarity cannot occur for large S. 

We conclude that multistationarity occurs in this motif. Note that statement (iii) 
implies that if multistationarity occurs for some values E j [i\ > F and some S,P, then 
increasing S eventually makes the system monostationary. Observe also that according to 
(i) and (ii) together, P m in is required to be so large that the condition F—Ej ii\ > P m in/o' f _ l 
is not fulfilled. 
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Motif (1). This motif is a combination of Motif (c) and Motif (j). The chemical reactions 
of the system are: 

Qi E E CL^ E 

S + E^X X ^S X + E P + S 1 ^X 2 ^P 1 + S 1 P + E^X 3 Xp 1 + E 

bf b$ bf 



zf C F , 



F 



S 1 +F 1 ^Y X ^S + Fx P X +F 2 ^Y 2 ^P + F 2 
bf q 

The differential equations describing the system are the following: 



Po 


= bfx 2 


+ bfX 3 


+ cfY 2 -afS 1 P -afEP 




So 


= bfX x + 


cfYi- 


afESo 




E 


= (bf+ 


cf)X x -\ 


- (bf + cf )X 3 - (afS + afP )E 




Xi 


= afES Q 


-(&f + 


cf)Xi 


(68) 




= (bf+ 


cf)Y, - 






x 2 


— afS\Po 


- (bf ~ 


f cf)X 2 


(69) 


F 2 


= 04+ 


cf)Y 2 - 


afF 2 P x 




x 3 


= afEPo 


+ 


ci)x 3 


(70) 


Pi 


= bfY 2 


+■ cfX 2 - 


hcfX 3 -afF 2 P x 


(66) 


Y x 


= afF x S x 


-(bf- 


Vcf)Y x 


(71) 


Si 


= cfX x 


+ bfY- 


f (bf + cf)X 2 - (afF x + afP )S x 


(67) 


Y 2 


= afF 2 P 1 


-ibf- 


Vc F 2 )Y 2 


(72) 



The conservation laws are given by 
E = E+X 1 +X 3 , F k = F k +Y k , S = S +S 1 +X 1 +X 2 +Y 1 , P = P 0+ P 1+ X 2 +Y 2 +X 3 . 
If total amounts are given, the system to be solved consists of 12 equations in 12 variables 



which are the equations for the total amounts and, for instance, equations (66)-(72). We 
obtain 

X 1 = r ]x ES< h X 2 = m S 1 P , X 3 = m EP , Y x = MiSi, Y 2 = 5 2 F 2 P X 
with constants rj k = af /(bf + cf) and 5k = af /(bf + cf). Equations (67), (69), (71 ), and 



(66), (72) provide 

X x = inY u ^X 2 + v 3 l X 3 - Y 2 = 

with \x x = cf /cf , fi 2 = cf /cf and ^l 3 = cf /cf . 

We write X x as a function of Y\ and isolate E,F X ,F 2 using the conservation laws for 
the total amounts E,F X ,F 2 . If E,Fk > then E, F& ^ at steady state and it follows 
that E — X x — X 3 > for any BMSS. Further, for any BMSS we require < Y k < F k and 
Y x < E / ji x . It follows as well that X 2 < fi 2 F 2 . If in addition, S > and P > then also 
S ,S x ,X 1 ,Y 1 ,X 2 ,P ,Pi,Y 2 ,X 3 > for any BMSS. 

Let £i = min(i ? i, E/fi X ), T x = (0, £i) and T 2 = (0, [i 2 F 2 ). It is convenient for this motif 
to exclude in the intervals. A necessary condition for positivity (i.e. a BMSS) is thus 
Y\ € Y X ,X 2 6 IV We proceed to write S x as an increasing C 1 functions of Y x £ T x and Pq 
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as a C 1 function increasing in X2 G r 2 and decreasing in Y\ G Ti: 

1 s 1 (F 1 -Y 1 y n mYi • 

Using X3 = t]^EPq, we express X3 as a positive C 1 function, increasing in X 2 G r 2 and 
decreasing in Yi G Ti: 

x 3 = ^(y 1 ,x 2) -^ Fl - y ^ E -^ 



772^ + 6 1 r j3 X 2 (F 1 - Y 1 ) 



We proceed to write Y2 as a positive C 1 function, increasing in X2 G r 2 and decreasing in 

Yi g r i: 

y 2 = $ y ( Yi , x 2 ) = iq 1 x 2 + ^3 1 $ x (Yi , x 2 ) . 

For Y 2 <F 2 , we require (Y 1 ,X 2 ) G r" with 

r' = {(Yi,x 2 )| Y t g r 1; x 2 g r 2 , $ y (Fi,x 2 ) < f 2 }. 

From Pi = y 2 /(5 2 -F 2 ) we have 

®y{Y 1 ,X 2 ) 



Pi = 



S 2 (F 2 -^ Y {Y 1 ,X 2 )) 



which is increasing in X 2 , decreasing in Y\ and positive and continuous provided (Yi, X2) G 
r'. Finally, 

s = Mi^i 

° m (E-^ 1 Y 1 -^ x (Y 1 ,X2)) 

which is positive and C 1 in T'. 

To sum up: All concentrations at steady state are positive if only if (Yi,X 2 ) G T' . To 
find a final relation between X 2 , Yi, we consider the total amount P. For Yi,X 2 in V, 

P = <p P (Yi,X 2 ), 

where (fp a positive C 1 function. Note also that for every fixed value of Y\ G T±, <pp(Yi, •) 
increases in X 2 , is well-defined at X2 = where it is zero and tends to infinity as <£y (Yi, X 2 ) 
tends to F 2 - Therefore, for P > and any Y\ G T\, there exists < X 2 satisfying 
P = tpp(Yi,X 2 ) and ^y(Y\,X 2 ) < F 2 . Thus, there exists a function / defined on Ti for 
which 

X2 = /(Yi) 

at steady state and (Yi,/(Yi)) G V. Since ipp is increasing in X2 and decreasing in Yi, 
the function / is C 1 and increasing in r := T\. 
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All concentrations at steady state are positive if and only if Y\ £ T. The concentrations 
S\ , X\ , Y\ are increasing functions of Y\ and independent of X2 ■ X2 is increasing in Y\ 
(after substitution of /). So might be increasing or decreasing in Y\ since &x(Yi, X2) is 
not increasing in Y\. In any case, if we insert these values into S, we obtain that the BMSSs 
are given by a relation 

S = V {Y 1 ), 

where cp is a positive C 1 function defined on T. When Y\ tends to £1, the function tends 
to infinity, since f(Y\) is bounded by n%Fi. When Y\ tends to zero, then X2 tends to zero 
as well as can be seen from the equations Y% = 5iFiS\ and X2 = ^Si-Po and the fact that 
Pq is bounded by P. Therefore, there is at least one BMSS. Multistationarity can only 
occur in the form of Figure [6^c), implying the existence of lower and upper bounds on S, 
Smin < S max , for the existence of multistationarity. 

If S*o is increasing in Y%, then so is tp and there is exactly one BMSS. The derivative of 
So with respect to Y\ is 

dSo = f i 1 (r f2 E + 6iV3X2(»iFi-E)) 8S , 

m mm (E - Ml y x ) 2 dX 2 J ' 

Since (dS /dX 2 )f' > 0, we have that if (a) ftFi - E > or (b) n x Fi - E < and 
r] 2 E + 5\ri^F2{^iFi — E) > 0, then §^ > and hence multistationarity cannot occur. 

Proposition 6. Fix the total amount F\. Then there exist total amounts P, S , F2, E 
satisfying fiiF\ — E < for which <f'(Yi) < for some Y\. 

It follows from the discussion above the proposition that 

7] 2 E + 61 7/3^2(^1^1 - E) < 

for multistationarity to occur. In fact, F2 and E are chosen so large that this condition is 
fulfilled. The proof of this proposition is provided in Appendix A. 

Result 12 (Motif (1)). Consider a cascade with different phosphatases acting in the two 
layers and where the kinase of the first layer also acts in the second layer. Assume that 
the total amounts E, F\, F2, S, P are positive. 

(i) If (a) n{Fi-E> or (b) m~F 1 -E< and r) 2 E + 5 1 r] 3 F2(viF 1 -E) > 0, then the 
system has exactly one BMSS. 

(ii) For any F\ and E,F2 large and satisfying fJ,\Fi — E < and r/2E + Sirj^F^ (mi-^i — 
E) < 0, there exist values P, S for which the system displays multistationarity. Fur- 
ther, in this case, there exists S m i n < S max for which multiple steady states occur for 
~R < S* < ~K 

" mm >-> *J max • 
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C Stability analysis 

We provide here the mathematical details of the stability analysis. We prove equation [3] 
in the main text and show that for the motifs exhibiting multistationarity, the equilibrium 
points for which ip' < are unstable. 



C.l The determinant of the Jacobian 

We prove here equation [3] of the main text. For any 1 < j < n, we define the projection 
ir^ of M n to R n_1 that removes the j-th coordinate by 

(xx , . . . , X n ) I—t- (xx , . . . , Xj , . . . , x n ) . 

For simplicity, let x^ = tt^\x) for any x € 1" and let Q^) = tt^\Q) C M n_1 for any open 
set Q C R n . Note that the latter is also an open set (projection maps are open maps). 

For a differentiable function / = (/i, . . . , /„) : C R n — y M n , we denote by J/(x) the 
Jacobian of / at x £ Q. The (z, j) entry of Jf(x) is dfi/dxj(x). 

In the next proposition, £1 is an open neighborhood of z, suitably chosen. 

Proposition 7. Let f = . . . ,/ n ): ^ C M n — >■ R n 6e a differentiable function defined 
on an open set Q, such that there is z 6 f2 u>i£/i /(z) = 0. Assume that Xj can be eliminated 
from the equation fi = 0; that is, there exists a differentiable function ip: C M ra_1 — > 
R such that Xj = ip(x^) on {x G Q\fi(x) = 0}. Dearie /: R n ~ l by f k (x) = 

fk{xi, . . . , Xj-i,ip(x),Xj, . . . , x n -i) for all k ^ i. Then, the determinant of the Jacobian of 
f at z satisfies 

t^ iz) det W zU) » = det(J/W). (73) 

Proof It is sufficient to show that the proposition holds for i = j = 1 . For other values of 
i and j the result is obtained by reorganizing rows and columns and keeping track of the 
sign of the determinant. For i = j = 1, the proposition states that 

y±{z) det(J f -(zW)) = det(J f (z)). 

where z^ = (z 2 , . . . , z n ) and z = (z\, ... ,z n ). 

If the Jacobian is written in column vector notation, we have that 

J f (z) = (D 1 f(z), D n f(z)) where D k f = . . . , ^ ' 



That is, D k f is the vector of the partial derivatives of / with respect to x k . 



51 



Enzyme sharing and multistationarity 



By assumption, z = {ip{z^), z%, . ..,z n ) G W l . Observe that if {dfi/dx\)(z) = then 
from ^(z)^(z^) + %±(z) = 0, we have that {dfi/dx k ){z) = for all k. Consequently, 
det(Jy(z)) = (the matrix has one row of zeroes) and hence the proposition holds. 

Hence, we can assume that {df\/dx\){z) / 0. By implicit differentiation, 

Additionally, by the chain rule we have {k ^ 1) 

M (z (i)) = ^L( Z) + W( Z) ?±( Z W\ 
dxk dxk dx\ dxk 

Let = (/2, . . . , / n ) be the projection of / onto the last n — 1 coordinates. From the 
equation above, it follows that 

Jj{z^) = (D 2 f, £>„/)(*«) = (D 2 fW, . . . , D n fW)(z) 



( 9 ^ /.im ^ raw .a 9 ^ 

where matrices are written as column vectors. 



Note that ^-(z^) is a scalar and Dif^\z)-§^-(z^) is the vector with Z-th component 



dx 

§^{z) ^^{z^). Further, using the multilinear expansion of a determinant, we obtain 
det(J f (z^)) = det((D 2 / (1) , . . . , D n f^)(z)) 

+ X>tCD 2 / (1) (*), • • • , D lf W(z)^% D n fW(z)) 



k=2 



+ ^2 det(a 2 ,...,a n ), 



{fci,...,fe ( }C{2,...,n} 

where a s = D s f( l \z) for s ^ k±, ... ,ki and a s = J^(-z^) otherwise. In each of 

these summands, there are at least two terms of the second form, say Dif^ (z)-^(z^) 
and Dif^ (z)-^^(z^) for some k / m. If the two columns are non-zero, they are linearly 
dependent and the determinant of each of the matrices (02, . . . , a n ) is zero. Further, note 
that 

det(Z? 2 /W(z), . . . , D lf U(z)^(zM), . . . , D n fW(z)) 

OXk 

= det(p 2 /«, . . . , D lf W,. . . , D n fW)(z)) 
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Finally, using (74), we obtain 



dfi 
dxi 



n 



dx k 



(z)det(J / (z( 1 )))=^(-l) fc - 



(z)det(D 1 f ( - 1 \...,D k fW,...,D n fW)(z) 



k=l 



= det(J/(z)) 



by considering the development of the determinant of Jf(z) along the first row. 



□ 



Application. Let a dynamical system in f2 C W 1 be given 



x = f{x) 



with x = (xi, . . . , x n ) and / = (/i, . . . , f n ). Let z = (zi, . . . , z n ) be an equilibrium point, 
i.e., f(z\, . . . , z n ) = 0. Let det(Jf(z)) be the determinant of the Jacobian of / at z. We 
make the following observation: 

• If n is odd and det(J/(z)) > 0, then z is unstable. 

• If n is even and det(J/(z)) < 0, then z is unstable. 

Indeed, if z is (asymptotically) stable and det(Jj(z)) / 0, all eigenvalues have negative 
real parts. Since det(J/(z)) is a real number (it is the determinant of a real matrix), 
the complex eigenvalues come in pairs of conjugates a + bi, a — bi and the product of the 
eigenvalues is a positive number. If n is odd and all eigenvalues have negative real parts, 
their product must be negative and hence the determinant of the Jacobian, which equals 
the product of the eigenvalues, must be negative. If n is even and z (asymptotically) stable, 
then the product of the eigenvalues must be positive. 

An equilibrium point satisfies f(zi,...,z n ) = 0. If elimination of variables can be 
applied then we can use Proposition [7] to track the sign of the determinant and potentially 
detect instability. 

Chemical reaction networks. Consider any of the motifs, or in general, any chemical 
reaction network with conservation laws. The conserved total amounts imply that the 
dynamics of the associated dynamical system takes place in a fixed subspace of W 1 . In 
general, we have a dynamical system 



x = f(x) 



and a series of (independent) conservation laws 



gi(x) - a = 0, 



i = 1 . . . , k 



(75) 
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where g% a linear function of x satisfying gi(x) = and Ci £ M. (these correspond to the 
total amounts). By independence we mean that the rank of this linear system is k. The 
conservation laws do not depend on the rate constants and gi does not depend on the total 
amounts. 

The existence of conservation laws implies that the determinant of the Jacobian of / 
at any point x is zero, since the matrix has linear relations among the rows. Therefore, 
stability of equilibrium points cannot be analyzed directly, but need to be considered inside 
the stoichiometry class they belong to. 

Since (75) is a linear system of rank k, Gauss elimination allows the elimination of k 
variables. For simplicity, we can rename the variables such that the eliminated variables 
are x\, . . . ,Xf, and those that remain are x k+ \, ■■■,x n . Apply the same renaming to the 
functions fi, such that if x\ is now variable Xj, function f\ is labeled fj. By elimination, 
there exist (polynomial) functions 

X-i — Tpi (Xfc+l , • • • , x n , c) , i — 1, . . . , k 

such that fi(ipi, . . . , ipk, Xk+i, • • • , x n ) = 0. Here c = (ci, . . . , Ck) is the vector of initial total 
amounts. 

For a fixed stoichiometry class c = (ci, . . . , Ck) and x = (xk+i, ■ ■ ■ , x n ), let 

fj(x,c) = fj(^i(x,c), . . .,ip k (x,c),x k+1 , . . .,x n ), j = k + 1, .. . ,n. 

To investigate stability of an equilibrium point (z\, . . . , z n ) belonging to the stoichiometry 
class c with Zi — ipi(zk+i, • • • , z n , c) = 0, we consider the the eigenvalues of the Jacobian of 
the function 

(fk+l{z,c), . . .,f n (z,c)), 

evaluated at z = (zk+i, ■ ■ ■ ,z n ). This function corresponds to the system called reduced in 
the main text. 

By Proposition [7] the sign of the determinant of this Jacobian at z is exactly the sign 
of the determinant of the Jacobian of the system 



Xi — ipi(x, c) = i — 1, . . . , k 
fi(x) = i = k + l,...,n 



evaluated at z. Indeed, the process leading from this system to the reduced system consists 
of successive eliminations with i = j = 1 and the derivative of the eliminated function 
corresponding to the eliminated variable is 1 (and thus positive). 

Note that the reduced system has n — k variables. Let J(Sf) denote the Jacobian of 
Sf and let z be an equilibrium point with total amounts c. We conclude that: 

Result 13. With the notation introduced above: 
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• If n — k is odd and det(J(Sf)(z)) > 0, then z is unstable. 

• If n — k is even and det(J(Sf)(z)) < 0, then z is unstable. 

C.2 Instability in the multistationary motifs 



Here we show how to apply Result 13 to each of the multistationary motifs: (f), (g), (i) 



(k), (1). Using Proposition [7J we find for all motifs but Motif (1) that 

sign((p' (zi)) = sign(det(J(5/)(z))), 
where z% is the variable of ip (typically an intermediate complex Y) and n — k is even. For 



Motif (1), n — k is odd and there is a change of sign. Hence, using Result 13, the steady 
state z is unstable whenever <p is decreasing. In particular, let e = ±1 (i.e. either e = +1 
or e = — 1) such that 

e • sign(^'(zi)) = sign(det(J(S/)(;s))). 

Remark 4. The sign of the determinant of a matrix remains unchanged if a linear com- 
bination of rows are added to another row (in fact the determinant does not change). If 
we multiply a row by a negative number, the determinant changes sign. For example, if 
the equation —(b E + c E )X + a E ESo = is transformed into X — rjESo = 0, then the sign 
of the determinant changes. If two such equations are transformed in this way, then the 
determinant remains with unchanged sign. In the sequel, the sign remains unchanged if the 
number of transformations of this type is even, or equivalently, if the number of constants 
5* , T]* is even. 

In the sequel, elimination is tracked using a table as in the main text and e = Jll^i 1 e i- 
The column 'Behavior' in the table shows (i,j, s) where i, respectively j, are the indices of 
the equation, respectively the variable, that iteratively are being eliminated and s indicates 
whether fi (/, after substitution of the previous eliminated variables) is increasing (s is +) 
or decreasing (s is — ) as function of Xj. 

Motif (f) is covered in the main text and it is thus skipped here. 

Motif (g). Consider the variables in the following order 

(xi, x 2 , x 3 , X4, x 5 , x G , x 7 , x$, x 9 ) = (E,F,S ,X 1 ,X 2 ,Si,S2,Y 2 ,Y 1 ). 

The variables of the reduced system are X\ , X 2 , S\ , S 2 , Y 2 , Y\ and the variables eliminated 
using the conservation laws are E,F,Sq. The sign of the determinant of the Jacobian of 
the system Sf is the same as the sign of the determinant of the Jacobian of the system 
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f l {x) = E + X l +X 2 -E 
f 2 (x) = F + Y 1 +Y 2 -F 

h(x) = S + S l + S 2 + X l +X 2 + Y l + Y 2 -S 
fi(x) = Xi - t]iESq 



f 5 (x) = X 2 -r ]2 ES 1 
h{x) = X X - mY x 
h{x) = X 2 - fi 2 Y 2 
h(x) = Y 2 -5 2 FS 2 
f 9 (x) = Y 1 -5 l FS 1 . 



With the notation introduced above, x = (Xi,X 2 , Si, S 2 , Y 2 , Yi) and 

ih(x) = -X 1 -X 2 +E, ^ 2 {x) = -Y l -Y 2 +F, ?p 3 (x) = -S 1 -S 2 -X 1 -X 2 -Y 1 -Y 2 + S. 
Here n = 9 and k = 3, such that n — k = 6 is even. The function (p is equation f% with 



all variables but xq = Si eliminated. The eliminations we performed in Section B.2| (Motif 
(g)) are summarized in the following table: 



I 


Elimination 


Behavior 




1 


Elimination 


Behavior 




1 


(h,E) 


(1,1,+) 


+ 


6 


(/4, So) 


(2,1,-) 


+ 


2 


(f2,F) 


(1,1,+) 


+ 


7 


(fs,S 2 ) 


(3,2,-) 


+ 


4 


(/6,*l) 


(4,2,+) 


+ 


8 


(/ 9 ,n) 


(3,3,+) 


+ 


5 


(fr,X 2 ) 


(4,2,+) 


+ 


9 


(h,Y 2 ) 


(2,2,+) 


+ 



Since e = 1, we conclude that sign((p' (zq)) = sign(det( J(Sf)(z))) for any equilibrium 
point z with zq = S\. 



Motif (i). Consider the variables in the following order 

(x 1 , x 2 , x 3 , x 4 , x 5 , x 6 , x 7 , x s , x 9 , x w ) = (E, F, So, Pq, X\,X 2 , S\,P\,Y\,Y 2 ). 

The variables of the reduced system are X\ , X 2 ,S\,P\,Y\,Y 2 and the variables eliminated 
using the conservation laws are E, F, So, Po- The sign of the determinant of the Jacobian 
of the system Sf is the same as the sign of the determinant of the Jacobian of the system 

f l {x) = E + X l +X 2 -E f 5 (x) = X 1 -r ]1 ESo f 8 (x)=X 2 -fi 2 Y 2 

f 2 {x) = F + Y l +Y 2 -F f & (x) = X 2 - m EP h (x) = Y 1 - StFSx 

f 3 (x) = So + S 1 + X 1 + Y 1 -S f 7 (x) = X 1 - f i 1 Y 1 fio{x)=Y 2 -5 2 FP l 
h(x) = Po + Pi+X 2 + Y 2 -P 
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Let x = (X 1 ,X 2 ,S 1 ,P 1 ,Y U Y 2 ). Then 
^i(x) = -X 1 -X 2 + E 
Mx) = -Y l -Y 2 + F 

4, such that n — k 



il> 3 (x) = -S 1 -X 1 -Y 1 + S 
Mx) = -Pi~X 2 -Y 2 + P. 

6 is even. The function <p is equation with all 



Here n = 10 and k 

variables but x±o = Y 2 eliminated. The eliminations we performed in Section B.3 (Motif 
(i)) are summarized in the following table: 



/ 


Elimination 


Behavior 




1 


Elimination 


Behavior e/ 


1 


(fl,E) 


(1,1,+) 


+ 


6 


(/io,Pi) 


(5,3,-) - 


2 


(f2,F) 


(1,1,+) 


+ 


7 


(h,s ) 


(3,1,-) 


3 


(fl,Xl) 


(5,3,+) 


+ 


8 




(3,1,-) - 


4 


(fs,X 2 ) 


(5,3,+) 


+ 


9 


(/3,H) 


(1,1,+) + 


5 




(5,3,-) 











Note that the last elimination f% corresponds to tpi(Yi, Y 2 ) — S = 0, which is increasing in 
Y\. Since e = Ilf=i e i = 1, we coi 
equilibrium point z and z\q = Y 2 . 



Y\. Since e = Ilf=i e i = 1> we conclude that sign(<//(,zio)) = sign(det(J(5/)(z))) for any 



Motif (k). Consider the variables in the following order 

(xi,x 2 ,x 3 ,X4,,XB,x e ,xr,xs,X9,xio) = (E, F, S , P , X 1 , X 2 , Si , Pi , Yi , Y 2 ) . 

The variables of the reduced system are X\ , X 2 , Si , Pi , Y\ , Y 2 and the variables eliminated 
using the conservation laws are E, F, So, Po. The sign of the determinant of the Jacobian 
of the system Sf is the same as the sign of the determinant of the Jacobian of the system 



f 1 {x) = E + X 1 -E 
f 2 (x) = F + Y l +Y 2 -F 
h(x) = S + S 1 + X 1 +Y 1 +X 2 
U(x) = P + Pi + X 2 + Y 2 - P 
Ux = (X 1 ,X 2 ,S 1 ,P 1 ,Y 1 ,Y 2 ) then 

^i(x) = -X l + E 
Mx) = -Yi-Y 2 + F 



S 



k(x) 
fr(x) 



Xi 

x 2 

X! 



ViESq 

mS\Po 



h{x) 

h{x) 

ho(x) 



X 2 - u 2 Y 2 
Y 1 - 6 1 FS 1 
Y 2 - 5 2 FP V 



Mx) = -S 1 -X 1 -Y 1 -X 2 + S 
Mx) = -Pi~X 2 -Y 2 + P. 

6 is even. The function (p is equation with all 



Here n = 10 and k = 4, such that n — k 
variables but x\q = Y 2 eliminated. The eliminations we performed in Section B.4 (Motif 
(k)) are summarized in the following table: 
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/ 


Elimination 


Behavior 




1 


Elimination 


Behavior e/ 


1 


ifl,E) 


(1,1,+) 


+ 


6 


(ho, Pi) 


(5,3,-) - 


2 


(f2,F) 


(1,1,+) 


+ 


7 


(h,s ) 


(3,1,-) - 


3 


(f7,Xl) 


(5,3,+) 


+ 


8 


(h,Po) 


(3,1,-) - 


4 


(f 8 ,x 2 ) 


(5,3,+) 


+ 


9 


(/ 3 ,n) 


(1,1,+) + 


5 


(f9,St) 


(5,3,-) 











Note that the last elimination f 3 corresponds to <&(Yx, Y 2 ) — S = 0. Since e = nf=i € i = 1, we 
conclude that sign(<£>' (z\o)) = sign(det( J(Sf)(z))) for any equilibrium point z and z±o = Y 2 . 

Motif (1). Consider the following order of the variables 

(x 1 ,x 2 ,x 3 ,X4,x 5 ,x G ,x 7 ,x s ,xc l ,x w ,xn,x 1 2) = (E,F 1 ,F 2 ,Sq,Pq,X 1 ,X 2 ,X 3 ,S 1 ,P 1 ,Y 1 ,Y 2 ). 

The variables of the reduced system are Xx, X 2 , X 3 , Si, Pi, Y±, Y 2 and the variables elimi- 
nated using the conservation laws are E, F\,F 2 ,Sq,Pq. The sign of the determinant of the 
Jacobian of the system Sf is opposite the sign of the determinant of the Jacobian of the 
system 



fi(x 
f±(x 

h{x 

fe(x 



E + X\ + ^3 — E 
Fx + Yx-F 
F 2 + Y 2 - F 2 

S + S 1 +X 1 +Y 1 +X 2 -S 
P + P 1 +X 2 + Y 2 + X 3 -P 
Xx - mESo 



f 7 (x) =X 2 - ri 2 SxP 
fs(x) =X 3 - m EP 
f 9 (x) = Xx-nxYx 
fxo(x) = X 2 - (fi 2 /fi 3 )X 3 - n 2 Y 2 

fxx(x) = Yx-8xFxSx 
f 12 (x)=Y 2 -S 2 F 2 Px. 



Indeed, the reason for the change in sign comes from Remark[4j since there is an odd number 
of equations that are rearranged by multiplication of negative numbers (corresponding to 

m^2,v3,Sx,s 2 ). 

Let x = (Xx,X 2 ,X 3 ,Sx,Px,Yx,Y 2 ). Then 



ipx (x) 
ip 2 (x) 
il) 3 {x) 



-Xx — X 3 + E 
-Yi + Fx 
-Y 2 + F 2 . 



ip 5 (x) 



-Sx 
Pi 



Xx-Yx-X 2 + S 
X 2 - Y 2 - X 3 + P 



Here n = 12 and k = 5, such that n — k = 7 is odd. In this case, we should see that e = 1. 
The function ip is equation with all variables but xxx = Yx eliminated. The eliminations 
we performed in Section B.4 (Motif (1)) are summarized in the following table: 
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/ 


Elimination 


Behavior 




1 


Elimination 


Behavior 




1 




(1,1,+) 


+ 


7 


(/8,*3) 


(4,3,+) 




2 


(f2,F 1 ) 


(1,1,+) 


+ 


8 


(ho,Y 2 ) 


(4,5,-) 


+ 


3 


(fs,F 2 ) 


(1,1,+) 


+ 


9 


Ul2,Pl) 


(4,3,-) 


+ 


4 


(fv,Xl) 


(6,3,+) 




10 


if 6, So) 


(3,1,-) 




5 


if n, Si) 


(7,5,-) 




11 


(h,x 2 ) 


(2,1,+) 




6 


(f 7 ,Po) 


(4,2,-) 













The last elimination /s corresponds to (fp(Yi, X 2 ) — 5 = 0, which is increasing in X 2 . 
Since e = Y\]Li e i = 1, we conclude that sign((^'(zio)) = — sign(det(J(5j)(z))) for any 
equilibrium point z and z\\ = Y\. Together with the fact that n — k is odd, the steady 
states for which ip' < are unstable. 

C.3 The monostationary motifs 

For Motifs (a)-(d), (e), (h) and (j), the above procedure is non-conclusive. We only show 
this for Motif (j), since for the other motifs we can prove that the steady state is asymp- 
totically stable using the Routh-Hurwitz criterion. 

Motif (j). Consider the variables in the following order 

(xi,x 2 ,x 3 ,x 4: ,x 5 ,xg,x 7 ,x 8 ,x 9 ,x 1 o,xi 1 ) = (E,F 1 ,F 2 ,S ,P ,Xi,X 2 ,Si,P 1 ,Y 1 ,Y 2 ). 

The variables of the reduced system are X\ , X 2 , Si , Pi , Y\ , Y 2 and the variables eliminated 
using the conservation laws are E, F\ , F 2 , So , Pq . The sign of the determinant of the Jaco- 
bian of the system Sf is the same as the sign of the determinant of the Jacobian of the 



system 


















= E + X l -E 


h(x) 


= Xi 


- m E So 


h{x) 


= x 2 


- V2Y2 


h{x) 


= F l +Y l -F l 


fv(x) 


= x 2 


- mSiPo 


fio(x) 


= Yi 


- S l F 1 S 1 


h{x) 


= F 2 + Y 2 -F 2 


fs(x) 


= Xi 


- Ml^l 


fn(x) 


= Y 2 


- 5 2 F 2 P! 


fdx) 


= S + S 1 +X 1 +Y 1 + X 2 -S 














h(x) 


= P + P 1 + X 2 + Y 2 -P 















Here n = 11 and k = 5, such that n — k = 6 is even. The function ip is equation f§ with all 
variables but xn = Y 2 eliminated. The eliminations we performed in Section B.4 (Motif 
(j)) are summarized in the following table: 
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I 


Elimination 


Behavior 




1 


Elimination 


Behavior 




1 


(h,E) 


(1,1,+) 


+ 


6 


(/6, So) 


(3,1,-) 




2 


(f2,F 1 ) 


(1,1,+) 


+ 


7 


(ho, Si) 


(4,2,-) 




3 


(h,F 2 ) 


(1,1,+) 


+ 


8 


(f 7 ,Po) 


(3,1,-) 




4 


(f*,Xi) 


(5,3,+) 


+ 


9 


(hi, Pi) 


(3,1,-) 




5 


(h,x 2 ) 


(5,3,+) 


+ 


10 


(h,Yi) 


(1,1,+) 


+ 



Since e = n«=i e i = 1, we conclude that sign(<p' (zn)) = sign(det( J(Sf)(z))) for any equi- 
librium point z and z\\ = Y 2 . Since n — k is even, stable steady states have positive 
determinant, and since <p is always increasing, nothing can be concluded. 



Routh-Hurwitz. We have computationally checked that Motifs (a)-(d), (e) and (h) have 
asymptotically stable BMSSs. The Routh-Hurwitz criterion establishes when all roots of a 
polynomial have real negative parts from a condition on the coefficients of the polynomial. 
Specifically, consider a real polynomial 



p(z) 



a z + a\z H V a n -iz + a r 



where we can assume that ckq > 0. One constructs the following matrix: 



A 



( oti a 3 a 5 

«o a 2 Oti 

OL\ Q!3 

ao OL2 



Q 6 
Q4 



o \ 





V 



That is, the entry (i,j) of the matrix is a 



a 



nj 



i+2(j-i) 



with the convention that if i + 2(j — i) > n 
or i + 2(j — i) < 0, then aj +2 (j-i) = 0- The criterion states that if all leading principal 
minors of the matrix have positive sign, then all roots of the polynomial p(z) have negative 
real parts. 

In our case, the polynomial to be analyzed is the characteristic polynomial of the 
Jacobian of the reduced system at a steady state. For example, by eliminating the variables 
E and So, the reduced system of Motif (b) is 



X = -(&! + Cl )X + ai (E - X - Y) (S 
Y = -(6a + c 2 )Y + a 2 (E -X- Y)Si 
Si = aX + b 2 Y - a 2 (E - X - Y)S 1 . 



Sr-X-Y) 
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The Jacobian can easily be computed using a program that handles symbolic compu- 
tations (e.g. Mathematica™). For Motif (b) the Jacobian J is 

-&! - cx - a x (E + S - Si - 2X - 2Y) -ai(E + S - Sj - 2X - 2Y) -a^E-X-Y) 
-a 2 Si -b 2 -c 2 -a 2 S\ a 2 {E-X-Y) 

ci + a 2 Si b 2 + a 2 Si —a 2 (E-X — Y) 

Note that for any BMSS, the terms E + S — S\ — 2X — 2Y and E — X — Y are positive. The 
leading principal minors of the Routh-Hurwitz matrix are polynomials in the entries of J. 
The task of computationally determining the sign of these minors is greatly simplified if 
we substitute back E = E + X + Y , and S = So + Sx+X + Y into J and write: 

-b\ — ci — cl\{E + Sq) —a\{E + So) 
J = | — (I2S1 —b 2 — C2 — a2<Si 

cx + a 2 Si b 2 + a 2 Sx 

The entries of the matrix J are now in the set of variables 

V = {E,X,Y, So,Sx,ax,a2,bx,b 2 ,cx,C2} 

and all variables in V take positive values at any BMSS. Thus, each of the entries of J has a 
fixed sign for any BMSS. The coefficients of the characteristic polynomial are polynomials 
in V and thus, the leading principal minors of the Routh-Hurwitz matrix are also polyno- 
mials in V. If the coefficients of these polynomials are all positive, then we are guaranteed 
that they take positive values for any choice of positive rates and any set of positive concen- 
trations. This computation can easily be implemented and checked with Mathematica™. 
If we used the matrix J involving the terms E + S — Si — 2X — 2Y and E — X — Y , then 
the corresponding polynomials take values in {E, S, X, Y, Sx, ax, a 2 , bx, b 2 ,cx, c 2 } but their 
coefficients are no longer positive. Positivity checking requires a convenient grouping of 
terms. 

This procedure shows stability for all motifs but Motif (j) in which case some negative 
terms in the expansion of some of the minors occur. These are computationally difficult to 
handle and we have not been able to show that the minors all are positive. However, by 
random generation of values, it seems that the steady state is stable, though it remains to 
be proven. 



D Proofs 

Proposition^ Figure [7] illustrates the different cases of the proposition. Recall that T = {Y\ € 
ri|/(Yi) e r 2 } and Ai = nxFx - E. The derivative of / is 

^ _ mjSxEFx - 28iiJ, 1 Fiy + (5x - m)Hiy 2 ) 
H 2 (S 1 (F 1 -y)+ri 2 y) 2 
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The denominator is always positive. The numerator is a degree two polynomial in y which might 
have positive real roots and is positive for y = 0. We consider y < £1 = min(i 7 ' 1 , E/fj,). 

The numerator evaluated at y = F\ is rj 2 {8\E — (8\ + r) 2 )[iiFi)Fi. It is positive if and only if 

A= (1 + W<$i)AtiFi ~E< 0. 

If this is the case, then £1 = F\ and the numerator is bounded from below by 772/^1(^1(^1 — Z/) 2 + 
V2(Fi — y 2 )) > 0. Hence both the numerator and f'(y) are positive for all y < £1. 

Therefore, for A < the function / is increasing for all yer 1 = [0, F x ) (see FigureJ^A)). since 
/(0) = 0, the condition /(Yi) G T 2 is equivalent to Y a < / _1 (6) and hence T = [0, min(Fi, / _1 (6))- 
This proves (i). 

To prove (ii), note that f'(Fi) < for A > 0. Further, f(E/m) < if and only if J5(<Ji - 772) - 
SifiiFi < 0. This is clearly the case when Ai > 0. IfAi<0, then using A > 0, it follows that 

E{5 1 - 772) - tfi/xiFi < r] 2 {-E + (nFi) < 0. 

In either case, f'(E/ni) < 0. It follows that if A > 0, then, whatever the sign of Ai, we have 
/'(£i) < and hence the numerator of /' has exactly one positive real root a in Ti. Thus, / is a 
function that increases up to y — a, and decreases after a (see Figure]?]). We have that f(E/m) = 
and hence, if Ai > 0, £1 = E//j,±, /(£i) = 0, and the function / starts and ends in zero in T±. 
Oppositely, if A x < then £1 = F\ and /(£i) = — A1//X2 > 0. These differences are depicted in 
Figure gB-C). 

The three cases (ii)(a)-(c) of the proposition follow from the determination of T. The idea is 
depicted in Figure[8] For /(Yi) G T2 we require that /(Yi) < £2- The maximal value of / in Ti is 
/(a) which is strictly smaller than E/fi 2 because [i\Y\ + [i 2 Y 2 — Mi^i + M2/(ii) = E. Thus, we 
have 

• If f{a) < F 2 , then for all Y x G Ti, we have /(Yi) < /(a) < £2 = min(F 2 , E2/H2). Thus, 
T = [0, £1), which corresponds to case (c) of the proposition. 

• If f{a) > F2, then the horizontal line Y2 = £ 2 = F 2 intersects the graph of / over [0,E/^,±) 
in two points corresponding to Yi-values a.\ < a < a 2 < E/fi\ (a\ = a 2 if £2 = /(«))■ We 
also have that ot\ < £1, since a\ < a < 

In the latter case, if Y\ e [0,ai) U (a 2 , Ej /ii), then f(Yi) e r 2 , while f(Yi) lies outside T 2 if 
Yi G [ai,a2]. Therefore, if f(a) > F 2 , we have 

r = ( [o, Ql )u(a23// J i))n[o,6) = {|°' ai ; iir I^ 1 ;" 2 

(J0,ai) U (a 2 ,£i) if £1 > «2- 

The first case happens if and only if £1 = F\ and F\ < a 2 . Since a < F\, this condition is 
equivalent to Ai < and — Ai//i 2 = f(F±) > f{a 2 ) — F 2 , or < E — (n\F\ + I12F2). This proves 
(a). 

The second case is equivalent to either £1 = E/fii or £1 = -Fi > a 2 - Proceeding as above, this 
is equivalent to < /-ti-Fi + ^F 2 — E, which proves case (b). □ 
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Proposition^ Assume that A > 0, that is Proposition [ljii) applies. Let (fi(Yi) = <^i(Yl, f{Yi)) so 
that <p(Yi) = (pxiYi) + ip 2 (f(Yi)). The derivative of <p with respect to Y± is 

^( y i) = l^( r i) + l^(/( 1 i))^( y i) Wlth l^ = ( 1 + ^) , 



The term ^p^- is the only one that depends on F 2 . Note in addition that this term is decreasing in 
F 2 for Y 2 <F 2 . 

The variables Si . Yi, X\ are all increasing functions of Y\ . As for So = — ;== — ^7-^ — , we 
have that the derivative of /iiYi + fi 2 f(Yi) is 

(mYi + 5 1 (F 1 -Y 1 )y 

which is positive for Y\ < Fi,E//ii. Thus, So is also increasing in Yi and hence g^-(Yi) > 0. 

By Proposition [I] (ii) (b) and the discussion following the proposition, if F 2 — f(a) there exist 
(many) values of Yi G T" for which (p'(Yi) < 0. Fix one such Yi. Since ip' is continuous in F 2 , then 
there exists an e > for which f'{Y\) < for any F 2 G (/(a) — e, f{a) + e). This ensures that for 
values of F 2 > /(a) close to f(a), we have multistationarity. 

On the other hand, the term §yK/(Yi)) i s positive for any Yi 6 I\ Therefore, if vanishes for 
some Yi, it must be that Yj > a, where / is decreasing. In the interval / = (a, £i), —df/dYi is pos- 
itive, bounded from above and independent of F 2 . It follows from the expression for df /dY± given 
in the proof of Proposition [T] Further, for Yi G /, dTp x jdYx is positive (stated after equation J76|)), 
bounded from below and independent of F 2 . The latter two statements follow from Result [ljand 



equation (76) upon differentiation (where ( |76| ) is used to differentiate So). 

Finally when F 2 tends to infinity, d(p 2 /dY 2 tends to zero uniformly for all Yi G /, since 
Y 2 = f(Yi) < f(a) is independent of F 2 . Putting it all together, we find that for F 2 large 

for all Yi G I. Consequently the function ip is increasing in all T and there cannot be multistation- 
arity. 

Therefore, for values of F 2 above f(a) but 'close' to f(a) there is multistationarity, and for 
large values of F 2 , there is monostationarity. If F 2 is such that ip(Yi) is increasing for some value 
Yi, then for any F 2 > F 2 , <p(Yi) must be increasing too (as ip'(Y"i) continues being positive when 



increasing F 2 , cf. (77l). 

This implies that the two regions (mono- versus multistationarity) are separated by a certain 
boundary value of F 2 , say M, where M is the infimum of all F 2 such that <p'(Yi) > for all 
Y\ G r. For F 2 = M, ip is also strictly increasing: The derivative ip' might be zero (as M is the 
infimum) but this can only happen in a finite number of points because the derivative is a non-zero 
rational function. Therefore, we see that for all F 2 > M, there is only one BMSS, while for all 
M > F 2 > /(a) there is multistationarity. □ 



Proposition^ Recall that a — (/Xi — ^2)(Mi^i 7 ?2 — ^2^2Vi) and ip — ((pi, <p 2 ). 
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Let us compute det( J^,(Y\, Y2)). We have, writing C F = F — Y1 — Y2 and Ce = E — fi\Y\ — [12X2, 
dpi F Y 2 fiiE 



dY 1 ^ s 1 cl 8iC% m c% 



dip 2 _ 1 , F Yx fi 2 E 



oy 2 ^ s 2 c F s 2 c% mcl 



E 



MlM2YJ 


dip! 




MlM2^1 


mc% 


dY 2 


SiC* 




iiiy.\Y x 


dip 2 


Y 2 


MlM2^2 




dY l 





Let C = — 1 — [L2 > and D = |^ — 1 — \i\ > 0. Then, the determinant is given as 

det(j^(yi,y 2 )) = (1 + mi)(i + M2) + (i + + (1 + ^2)^ + — ^-3 + - + 4 + s 



with 

ijq(ef- fjbiYiF - Y2E) M2 /, , y , M 2 y 2 , (M2-Mi)yy 2 



A = 



hmClCp S im C E C F \ C F C e C e C_ 



>F 



) 



B = ih{e f - ^f - y x e) _ mi (, , y 2 , Miy , (Mi-M2)yy 2 ^ 



6 2 ViC%C F 6 2 ViCeC f 



+ + Mj^i + (Mi - M2)r 1 r 2 \ 
V C^? Ce CeCf ' 



where in the last equality of both equations, we substitute F by Cp + Yi + Y 2 and E by Ce + 
M1Y1 + M2Yz- Hence, det( J^{Y\, Y 2 )) is a sum of positive terms together with 



yy 2 

SiS 2 ViV2C E C F ' 



N = X X Z - r<2 n 2 (M2^2^l(M2 - Mi) + Mi^i%(Mi - M2)) 



Therefore, if 

(T = (pi - M2)(Ml^l ? ?2 - M2<W) > 

then det( J^,(Y\, Y 2 )) > and the function ip is always increasing. This proves (i). 

Let us assume now that a < 0. We make the following observation: the only negative term 
N of det( J^(Yi, Y2)) is also the term with denominator of highest degree in Ce,Cf- Further, 
all numerators in the expression of det(J^(Yi,Y 2 )) above are bounded. Thus, by letting Ce,Cf 
simultaneously tend to zero, we can obtain a negative determinant. 

Since Cf = F — Y\ — Y 2 and Ce = E — [i\Y\ — H2Y2, it is possible to make them small 
simultaneously by varying Yi,Y 2 , if and only if the two lines rp\ F — Yi — Y 2 = and r# : E — 
Miy — M2y2 = intersect for valid values of Yi,Y 2 . The intersection of the two lines is the point 

(y,y 2 ) with 

m F-E ^F-E 

i\ — , 12 — ■ 

M2 - Mi Mi - M2 

If Mi - M2 > 0, then for positive intersection values Yi, Y 2 we require ^F < E < n\F. If ^ 2 — Mi > 0) 
then we require H\F < E < [I2F . 

Hence, if these conditions above are satisfied, we are guaranteed the existence of values of Yi, Y2 
for which the determinant is negative and thus multistationarity occurs. Further, if Ce,Cf are 
small then S, P must be large. 

Assume finally that a < but also one of the two following conditions holds: 
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(a) hi — fj,2 > and either fi 2 F > E or E > HiF. 

(b) fj,2 — Mi > an d either /xiF > E ov E > \i 2 F . 

In this case multistationarity cannot occur. Assume that (a) holds (case (b) is similar). Since 
Ml — M2 > 0, we have \iiF > ju, 2 F and: 

• If /iiF > fj, 2 F > E, then using —Y^E > — /UfeYfc-F, we have 

EF - fiiYiF - Y 2 E > ~EF - (J,iYF - fi 2 Y 2 F = FC E > 

and 

EF - YiE - H2Y2F >EF- ^YiF - \x-2?%F = FC E > 0. 

• If E > fiiF > H2F, then using —fikYkF > —Y^E, we have 

EF - nxY{F - Y 2 E >EF- YiE - Y 2 E = EC F > 

and 

EF - YiE - \x 2 Y 2 F >KF~ Y{E - Y 2 E = EC F > 0. 
Therefore, both A and B are positive and the proposition is proved. □ 



Proposition [JJ Recall that / is a function of Y 2 defined in equation ( 65 ) by elimination of Yi from 
S = $(Yi,Y" 2 )- Thus, we have that /' = -$ 2 /$i, where 

$2 := d§/dY 2 = ri 2 + — _ Yl > 0, 

6 1 (F-Y 1 -Y 2 y 

$1 := d<P/dYi = 1 + ui + — =J=± + — = > 0. 

' * m {E-HiYiY 5i{F-Yi-Y 2 y 



Recall that 



dP ^2<5i / f ov \ e i v v nn , dPi F - f + f'Y 2 

;{f{F - / -2Y 2 ) - / Y 2 (F -Y 2 )), and 



dY 2 m p yjy J " J ^ dY 2 S 2 (F-f-Y 2 ] 



2 



§f| < only if F - f + f'Y 2 < 0, that is 

— — ui E(F — Yi) F 

> $r(F - Yi) - <£ 2 Y 2 = (1 + - Yi) - m Y 2 + + — (78) 

- mi^i) di(F-Yi-Y 2 ) 

Since the last two summands are positive for valid values oiYi,Y 2 , a necessary condition for < 
is CI : (1 + Hi)(F — Yi) — fJ, 2 Y 2 < 0. Note that if 1 + fii > fi 2 , then the condition cannot be fulfilled 
and Pi is always an increasing function of Y 2 for any E,F,S. 

Since /' < 0, we require C2 : F — Yi — 2Y 2 < in order to have ||| < 0. 

It follows that C1,C2 are necessary conditions for multistationarity. If P,S are given, the 
following are restrictions on Yi,Y 2 : 

Yi+Y 2 <F, fiiYi<E, {l + m)Yi+H 2 Y 2 <S, \x 2 Y 2 < P. 
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Let o = min(S l , P). We find that fi 2 Y 2 < a and hence 

F -Y 1 -2Y 2 >F - E/fn - 2o/ f i 2 and (1 + m)(F - Y t ) - fi 2 Y 2 >{l + m)(F- E/m) - a. 

Let Op = min(l + £ii,/x 2 /2). H F > E/n± and o^(F — E/fi±) > o, then we are guaranteed that 
both C1,C2 > and multistationarity cannot occur. 

□ 

Proposition^ Note that for <J> (i.e. S) to tend to +oo for any fixed Y 2 ET 2 = [0, mm(S//i 2 , F)) C 
[0, F), either Yi must tend to Ej[i\ or to F — Y 2 . The first case arises if and only ifY 2 < F — E/fii. 

(i) Assume that Y 2 > F — E//j,i. Consider the functions $ 2) 3>i from the proof of Proposition [4] 
and the expression for |^ and The derivative /' = —$2/^1 tends to —1 as S tends to infinity 

(and Yi tends to F - Y 2 ). Likewise, it follows that f(F - f - 2Y 2 ) - f'Y 2 (F - Y 2 ) tends to zero, 



and hence Sfl 1 also tends to zero. From equation (78) in the previous proof it follows that jjrr > 



as Y\ approaches F — Y 2 , that is as S becomes large (the last term is positive and dominates the 

dPp 1 3Pi 
dY 2 f dY 2 



other terms). Hence the derivative of tp{Y 2 ) becomes positive: (p'(Y 2 ) = 1 + fi 2 + SIt 1 + |& > 0. 



This proves (i). 

(ii) Assume that Ai = niF — E > and consider Y 2 < F — E/fii. As S tends to infinity, Y\ 
tends to E/fii. The limit curve ipoo is 

<Poo{Y 2 ) = (1 + 1*2)12 + 



^E/^ tf 2 (Ai-y 2 ) 

for Ai = A1//X1. Existence of values of Y 2 for which ip'^iY^ < is equivalent to 

p(Y 2 ) :=<J 2 (Ai -Y 2 ) 2 (1 + M2 + /3(A 1 -2y 2 ))+Ai <0 

with /3 = fiifi 2 Si/ri 2 E. The function p(y 2 ) is a polynomial of degree 3. For Y 2 = 0, it is positive 
and when y2 tends to infinity, the polynomial tends to —00. Therefore, it takes negative values in 
Y 2 E [0, Ai) if and only if there is a root in this interval. At Y 2 = Ai it is positive, and hence, there 
is at least one root after Ai. The derivative with respect to Y 2 is 

P '{Y 2 ) = -2<5 2 (A 1 - y 2 )(l + fi 2 + /?(2Ax - 3y 2 )). 

One zero of the derivative is Y 2 = A 1; while the other is 7 = (1 + fi 2 + 2/3Ai)/3/3. Note that at A 1; 
p(Ai) is positive. Therefore, for negative values of p to occur between and Ai, we need 7 < Ai 
and p(j) < 0. 

We have 7 < Ai if and only if Ai > (1 + fi>2)/P, in which case Ai is a maximum and 7 is a 
minimum. Evaluating p in 7 gives the following condition: 

A 2 = 27^ lf i 2 2 S 2 1 7 l2 E{fi 1 F -E)- hiSuMWxF - (Mi + t? 2 (1 + n 2 ))Ef < 0, 

which concludes the proof of (ii). 

□ 
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Proposition^ It is convenient to write ip as a function of X%,Yx, i.e., without substituting X 2 by 
f(Yi). The derivative of (p with respect to Y\ takes the form 

, ^i(mE + SimX2{^Fi-E)) Fx / Miiftgr^) , A 

</? = = h 1 + /ii H = V — hi/ 

with only the first term susceptible of being negative since / is an increasing function. If we can 
control /' by varying some total amounts we might be able to make the derivative tp' negative. 

We will see that there exist values of X2, F 2 , and E for which cp' as a function of the two variables 
Y\ and X2 is negative. Then, by defining P = ipp(Yi, X 2 ), we obtain regions of multistationarity. 
Note that this is equivalent to the approach taken in Motif (i), when we studied the Jacobian of 
(<ps, fp) rather than the derivative of ip. 

Using the implicit function theorem, the derivative of / with respect to Y\ is /' = — §^j§x\ ■ 
We compute the partial derivatives of Pq, Pi, X 2 ,Y 2 , X3 for each term in the numerator and the 
denominator. First of all we repeat some definitions, 

X S> (V Y \ hriMFi-Y^E-mYj 



and 



It follows that 



V2Yi+5 m X 2 (F 1 -Y 1 ) ' 
Y 2 = * Y (Y 1 ,X 2 ) = ^ 2 1 X 2 + ^ 1 ^ X - 



dcpp _ djPp + Px) -1,,, , -i^x d<p P _ d(P + P 1 ) , n , _ x d$ x 

dX 2 ~ 8X 2 + ^ + ^ + ^3 )qx 2 , m - g Yi +U+M 3 )~gyr-. 

For the denominator we have: 

dP = ftgV-Yi) d$ x = ^773772^1(^1 - Y^JE - ^Y,) 

dX 2 772^! ' dX 2 ( % Yi +S m X 2 (F 1 -Y 1 )) 2 ' 
dPi _ F 2 

9X 2 " <5 2 (F 2 -$ y ) 2 ' lM2 + A * 3 dX 2 y 
For the numerator, we have 

dP SxFxXi d<$> x Sm^miFi -Y t ) V2Fi(E - /xiFi) 



0Ti r^ 2 9Y a 7 72 y 1 +5 1 7 ?3 X 2 (F 1 -y 1 ) (7? 2 Yi + Si^CFx - Fx)) 2 ' 

9Pi _ -F 2 d<f> x 

~dY ~ d 2/ x 3 (F 2 -$y) 2 W 

Now fix Yj, Fx and let X 2 , F 2 ,E tend to +00. Assume that F goes slower than X 2 , for example 
put X 2 — E for some p > 0. Let F 2 — $y = K be constant such that F 2 = K + /j 2 X 2 + ^ 3 $x 
and X 2 go at the same rate towards infinity. Note that we can vary X 2 ,F 2l E without changing 
Yi , Fx or restricting their range of definition. 
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With these assumptions the numerator of /', — ^r-, goes to infinity proportionally to X 2 . Also 
the denominator of /', goes to infinity proportionally to X 2 . Consequently, /' converges 

towards a constant that depends on F\ , Y\ and K (in addition to the rate constants) 

J ~^ K T72 1 • 

V2Y? M3 

Going back to the expression for ip' , we find that under the same assumptions, 

, -Hi6 1 r} 3 X 2 Fi , 

where terms that eventually vanish are not shown. It follows that for large E, X 2 = E and 
F 2 = K + n 2 1 X 2 + ^ 3 " 1 $x, the derivative ip' eventually becomes negative for any choice of Y\ and 
F\, as desired. 

□ 



